回报时间序列的预期回报和协方差

Question

我正在尝试模拟此处定义的 Matlab ewstats 函数：

https://it.mathworks.com/help/finance/ewstats.html

matlab给出的结果如下：

> ExpReturn = 1×2
0.1995    0.1002

> ExpCovariance = 2×2
0.0032   -0.0017
-0.0017    0.0010

我正在尝试使用 RiskPortfolios R 包复制该示例：

https://cran.r-project.org/web/packages/RiskPortfolios/RiskPortfolios.pdf

我使用的 R 代码是这个：

library(RiskPortfolios)

rets <- as.matrix(cbind(c(0.24, 0.15, 0.27, 0.14), c(0.08, 0.13, 0.06, 0.13)))
w <- 0.98

rets
w
meanEstimation(rets, control = list(type = 'ewma', lambda = w))
covEstimation(rets, control = list(type = 'ewma', lambda = w))

均值估计与示例中的相同，但协方差矩阵不同：

> rets
     [,1] [,2]
[1,] 0.24 0.08
[2,] 0.15 0.13
[3,] 0.27 0.06
[4,] 0.14 0.13
> w
[1] 0.98
> 
> meanEstimation(rets, control = list(type = 'ewma', lambda = w))
[1] 0.1995434 0.1002031
> 
> covEstimation(rets, control = list(type = 'ewma', lambda = w))
             [,1]         [,2]
[1,]  0.007045044 -0.003857217
[2,] -0.003857217  0.002123827

我错过了什么吗？ 谢谢

Answer 1

他们使用不同的算法。来自 RiskPortfolio 手册：

ewma ... 参见 RiskMetrics (1996)

来自 Matlab 帮助页面：

ewstats 函数与用于确定回报时间序列的预期回报和协方差的 RiskMetrics® 方法之间没有关系。

很遗憾，Matlab 没有告诉我们使用的是哪种算法。

Answer 2

如果使用type = "lw"，他们会给出相同的答案：

round(covEstimation(rets, control = list(type = 'lw')), 4)

##   0.0032 -0.0017
##  -0.0017  0.0010

Answer 3

对于那些最终需要在 R 中等效的 ewstats 函数的人，这里是我编写的代码：

ewstats <- function(RetSeries, DecayFactor=NULL, WindowLength=NULL){

  #EWSTATS Expected return and covariance from return time series.  
  #   Optional exponential weighting emphasizes more recent data.
  # 
  #   [ExpReturn, ExpCovariance, NumEffObs] = ewstats(RetSeries, ...
  #                                           DecayFactor, WindowLength)
  #  
  #   Inputs:
  #     RetSeries : NUMOBS by NASSETS matrix of equally spaced incremental 
  #     return observations.  The first row is the oldest observation, and the
  #     last row is the most recent.
  #
  #     DecayFactor : Controls how much less each observation is weighted than its
  #     successor.  The k'th observation back in time has weight DecayFactor^k.
  #     DecayFactor must lie in the range: 0 < DecayFactor <= 1. 
  #     The default is DecayFactor = 1, which is the equally weighted linear
  #     moving average Model (BIS).  
  #
  #     WindowLength: The number of recent observations used in
  #     the computation.  The default is all NUMOBS observations.
  #
  #   Outputs:
  #     ExpReturn : 1 by NASSETS estimated expected returns.
  # 
  #     ExpCovariance : NASSETS by NASSETS estimated covariance matrix.  
  #
  #     NumEffObs: The number of effective observations is given by the formula:
  #     NumEffObs = (1-DecayFactor^WindowLength)/(1-DecayFactor).  Smaller
  #     DecayFactors or WindowLengths emphasize recent data more strongly, but
  #     use less of the available data set.
  #
  #   The standard deviations of the asset return processes are given by:
  #   STDVec = sqrt(diag(ECov)).  The correlation matrix is :
  #   CorrMat = VarMat./( STDVec*STDVec' )
  #
  #   See also MEAN, COV, COV2CORR.


  NumObs <- dim(RetSeries)[1]
  NumSeries <- dim(RetSeries)[2]

  # size the series and the window
  if (is.null(WindowLength)) {
    WindowLength <- NumObs
  }

  if (is.null(DecayFactor)) {
    DecayFactor = 1
  }


  if (DecayFactor <= 0 | DecayFactor > 1) {
    stop('Must have 0< decay factor <= 1.')
  }


  if (WindowLength > NumObs){
    stop(sprintf('Window Length #d must be <= number of observations #d',
                 WindowLength, NumObs))
  }


  # ------------------------------------------------------------------------
  # size the data to the window
  RetSeries <- RetSeries[NumObs-WindowLength+1:NumObs, ]

  # Calculate decay coefficients
  DecayPowers <- seq(WindowLength-1, 0, by = -1)
  VarWts <- sqrt(DecayFactor)^DecayPowers
  RetWts <- (DecayFactor)^DecayPowers

  NEff = sum(RetWts) # number of equivalent values in computation

  # Compute the exponentially weighted mean return
  WtSeries <- matrix(rep(RetWts, times = NumSeries),
                     nrow = length(RetWts), ncol = NumSeries) * RetSeries

  ERet <- colSums(WtSeries)/NEff;

  # Subtract the weighted mean from the original Series
  CenteredSeries <- RetSeries - matrix(rep(ERet, each = WindowLength),
                                       nrow = WindowLength, ncol = length(ERet))

  # Compute the weighted variance
  WtSeries <- matrix(rep(VarWts, times = NumSeries),
                     nrow = length(VarWts), ncol = NumSeries) * CenteredSeries

  ECov <- t(WtSeries) %*% WtSeries / NEff

  list(ExpReturn = ERet, ExpCovariance = ECov, NumEffObs = NEff)
}