Demo entry 6637246

R

   

Submitted by anonymous on Aug 30, 2017 at 16:02
Language: S. Code size: 16.0 kB.

#-----------------第二章function VAR()-------------------
function (x, p = 1, output = T, include.mean = T, fixed = NULL) 
{
  #p是VAR(p),意思是t时刻的数据与它前面p个时刻数据有关
  if (!is.matrix(x)) 
    x = as.matrix(x)
  Tn = dim(x)[1]    #---有多少行
  k = dim(x)[2]     #---有多少列,列数就是k表示k维数据
  if (p < 1)        
    p = 1
  idm = k * p
  ne = Tn - p       #一共有ne个数据可以计算(预测)
  ist = p + 1       #从ist个数据开始是可以预测的,因为ist前面没有p个数据
  y = x[ist:Tn, ]
  # 下面到27行,将元数据组织起来得到书上的X,每一行是[x(t),x(t-1)...x(t-p)]
  if (include.mean) {
    idm = idm + 1
    xmtx = cbind(rep(1, ne), x[p:(Tn - 1), ])
  }
  else {
    xmtx = x[p:(Tn - 1), ]
  }
  if (p > 1) {
    for (i in 2:p) {
      xmtx = cbind(xmtx, x[(ist - i):(Tn - i), ])   #一共是p*k列
    }
  }
  
  ndim = ncol(xmtx)                        #ndim表示xmtx有多少列
  if (length(fixed) == 0) {
    paridx = matrix(1, ndim, k)
  }
  else {
    paridx = fixed
  }
  res = NULL
  beta = matrix(0, ndim, k)
  sdbeta = matrix(0, ndim, k)
  npar = 0
  for (i in 1:k) {
    idx = c(1:ndim)[paridx[, i] == 1]       #返回的是paridx中第i列中数字是1数据的行敿
                                            #c(1:10)[c(11:20)>15] 的返回值是6 7 8 9 10
    resi = y[, i]
    if (length(idx) > 0) {
      xm = as.matrix(xmtx[, idx])           #【可能? ## 默认下xm就是xmtx
      npar = npar + dim(xm)[2]
      xpx = t(xm) %*% xm                    #xpx是一个(p*k) * (p*k)矩阵?
      xpxinv = solve(xpx)                   #求矩阵的逆
      xpy = t(xm) %*% as.matrix(y[, i], ne, 1)  #xpy是一个p*k * 1矩阵     【注】p*k是xmtx的列数 = ndim
      betai = xpxinv %*% xpy                #betai= (t(xm) . xm)^(-1) . t(xm) . t(y[,i]) 
      beta[idx, i] = betai                  #betai是一个p*k * 1
      resi = y[, i] - xm %*% betai          #resi ne*1
      #以上用到公式:P49定理2.2 vec(beta^) 所以说上面方法用的是ML估计 
      nee = dim(xm)[2]                      # nee = npar = ne
      sse = sum(resi * resi)/(Tn - p - nee) # 下面几行与模型检验有关,这行揭示残差协方差矩阵与残差(residual)之间的关系
      dd = diag(xpxinv)
      sdbeta[idx, i] = sqrt(dd * sse)
    }
    res = cbind(res, resi)
  }
  sse = t(res) %*% res/(Tn - p)             # 公式 2-55 求残差协方差矩阵
  aic = 0
  bic = 0
  hq = 0
  Phi = NULL
  Ph0 = NULL
  jst = 0
  if (include.mean) {
    Ph0 = beta[1, ]
    se = sdbeta[1, ]
    if (output) {
      cat("Constant term:", "\n")
      cat("Estimates: ", Ph0, "\n")
      cat("Std.Error: ", se, "\n")
    }
    jst = 1
  }
  if (include.mean) {
    for (i in 1:k) {
      if (abs(Ph0[i]) > 1e-08) 
        npar = npar - 1
    }
  }
  if (output) 
    cat("AR coefficient matrix", "\n")
  for (i in 1:p) {
    phi = t(beta[(jst + 1):(jst + k), ])
    se = t(sdbeta[(jst + 1):(jst + k), ])
    if (output) {
      cat("AR(", i, ")-matrix", "\n")
      print(phi, digits = 3)
      cat("standard error", "\n")
      print(se, digits = 3)
    }
    jst = jst + k
    Phi = cbind(Phi, phi)
  }
  if (output) {
    cat(" ", "\n")
    cat("Residuals cov-mtx:", "\n")
    print(sse)
    cat(" ", "\n")
  }
  dd = det(sse)
  d1 = log(dd)
  aic = d1 + (2 * npar)/Tn
  bic = d1 + log(Tn) * npar/Tn
  hq = d1 + 2 * log(log(Tn)) * npar/Tn
  if (output) {
    cat("det(SSE) = ", dd, "\n")
    cat("AIC = ", aic, "\n")
    cat("BIC = ", bic, "\n")
    cat("HQ  = ", hq, "\n")
  }
  VAR <- list(data = x, cnst = include.mean, order = p, coef = beta, 
              aic = aic, bic = bic, hq = hq, residuals = res, secoef = sdbeta, 
              Sigma = sse, Phi = Phi, Ph0 = Ph0)
}


#-----------------------------------VARpre函数---------------------------------------------------
function (model, h = 1, orig = 0, Out.level = F) 
{
  x = model$data
  Phi = model$Phi
  sig = model$Sigma    #Sigma是一个正定协方差矩阵,参考定理2.2里面的定义
  Ph0 = model$Ph0  
  p = model$order      #var(p)
  cnst = model$cnst    #ture
  np = dim(Phi)[2]
  k = dim(x)[2]
  nT = dim(x)[1]
  k = dim(x)[2]
  if (orig <= 0) 
    orig = nT
  if (orig > nT) 
    orig = nT
  psi = VARpsi(Phi, h)$psi 
  beta = t(Phi)
  if (length(Ph0) < 1) 
    Ph0 = rep(0, k)
  if (p > orig) {
    cat("Too few data points to produce forecasts", "\n")
  }
  pred = NULL
  se = NULL
  MSE = NULL
  mse = NULL
  px = as.matrix(x[1:orig, ])  #orig用于数据分割,只是取前orig行。
  Past = px[orig, ]
  if (p > 1) {
    for (j in 1:(p - 1)) {
      Past = c(Past, px[(orig - j), ])   #past是px的后p行展开成一维向量 即(px[-1],px[-2],...,px[-p])
    }
  }
  #---看chapter.R上有输出和一些解释
  cat("orig ", orig, "\n")
  ne = orig - p
  xmtx = NULL
  P = NULL
  if (cnst) 
    xmtx = rep(1, ne)
  xmtx = cbind(xmtx, x[p:(orig - 1), ])
  ist = p + 1
  if (p > 1) {
    for (j in 2:p) {
      xmtx = cbind(xmtx, x[(ist - j):(orig - j), ])
    }
  }
  xmtx = as.matrix(xmtx)
  G = t(xmtx) %*% xmtx/ne
  Ginv = solve(G)
  P = Phi
  vv = Ph0
  if (p > 1) {
    II = diag(rep(1, k * (p - 1)))
    II = cbind(II, matrix(0, (p - 1) * k, k))
    P = rbind(P, II)
    vv = c(vv, rep(0, (p - 1) * k))
  }
  if (cnst) {
    c1 = c(1, rep(0, np))
    P = cbind(vv, P)
    P = rbind(c1, P)
  }
  Sig = sig
  n1 = dim(P)[2]
  MSE = (n1/orig) * sig
  for (j in 1:h) {
    tmp = Ph0 + matrix(Past, 1, np) %*% beta
    px = rbind(px, tmp)
    if (np > k) {
      Past = c(tmp, Past[1:(np - k)])
    }
    else {
      Past = tmp
    }
    if (j > 1) {
      idx = (j - 1) * k
      wk = psi[, (idx + 1):(idx + k)]
      Sig = Sig + wk %*% sig %*% t(wk)
    }
    if (j > 1) {
      for (ii in 0:(j - 1)) {
        psii = diag(rep(1, k))
        if (ii > 0) {
          idx = ii * k
          psii = psi[, (idx + 1):(idx + k)]
        }
        P1 = P^(j - 1 - ii) %*% Ginv
        for (jj in 0:(j - 1)) {
          psij = diag(rep(1, k))
          if (jj > 0) {
            jdx = jj * k
            psij = psi[, (jdx + 1):(jdx + k)]
          }
          P2 = P^(j - 1 - jj) %*% G
          k1 = sum(diag(P1 %*% P2))
          MSE = (k1/orig) * psii %*% sig %*% t(psij)
        }
      }
    }
    se = rbind(se, sqrt(diag(Sig)))
    if (Out.level) {
      cat("Covariance matrix of forecast errors at horizon: ", 
          j, "\n")
      print(Sig)
      cat("Omega matrix at horizon: ", j, "\n")
      print(MSE)
    }
    MSE = MSE + Sig
    mse = rbind(mse, sqrt(diag(MSE)))
  }
  cat("Forecasts at origin: ", orig, "\n")
  print(px[(orig + 1):(orig + h), ], digits = 4)
  cat("Standard Errors of predictions: ", "\n")
  print(se[1:h, ], digits = 4)
  pred = px[(orig + 1):(orig + h), ]
  cat("Root mean square errors of predictions: ", "\n")
  print(mse[1:h, ], digits = 4)
  if (orig < nT) {
    cat("Observations, predicted values,     errors, and MSE", 
        "\n")
    tmp = NULL
    jend = min(nT, (orig + h))
    for (t in (orig + 1):jend) {
      case = c(t, x[t, ], px[t, ], x[t, ] - px[t, ])
      tmp = rbind(tmp, case)
    }
    colnames(tmp) <- c("time", rep("obs", k), rep("fcst", 
                                                  k), rep("err", k))
    idx = c(1)
    for (j in 1:k) {
      idx = c(idx, c(0, 1, 2) * k + j + 1)
    }
    tmp = tmp[, idx]
    print(round(tmp, 4))
  }
  VARpred <- list(pred = pred, se.err = se, mse = mse)
}


#--------------VMA function------------
function (da, q = 1, include.mean = T, fixed = NULL, beta = NULL, 
          sebeta = NULL, prelim = F, details = F, thres = 2) 
{
  if (!is.matrix(da)) 
    da = as.matrix(da)
  nT = dim(da)[1]
  k = dim(da)[2]
  if (q < 1) 
    q = 1
  kq = k * q
  THini <- function(y, x, q, include.mean) {
    if (!is.matrix(y)) 
      y = as.matrix(y)
    if (!is.matrix(x)) 
      x = as.matrix(x)
    nT = dim(y)[1]   # 610
    k = dim(y)[2]    # 2
    ist = 1 + q
    ne = nT - q
    if (include.mean) {
      xmtx = matrix(1, ne, 1)
    }
    else {
      xmtx = NULL
    }
    ymtx = y[ist:nT, ]
    for (j in 1:q) {
      xmtx = cbind(xmtx, x[(ist - j):(nT - j), ])
    }                                                               ## xmtx: ne*(kq+1)
                                                                    ## ymtx: ne*k
    xtx = crossprod(xmtx, xmtx)  #crossprod(x,y)相当于t(x) %*% y    ## xtx : (kq+1)*(kq+1)
    xty = crossprod(xmtx, ymtx)                                     ## xty : (kq+1)*k
    xtxinv = solve(xtx)         #xtx的逆
    beta = xtxinv %*% xty                                           ## beta: (kq+1)*k
    resi = ymtx - xmtx %*% beta 
    sse = crossprod(resi, resi)/ne
    dd = diag(xtxinv)
    sebeta = NULL
    for (j in 1:k) {
      se = sqrt(dd * sse[j, j])
      sebeta = cbind(sebeta, se)
    }
    THini <- list(estimates = beta, se = sebeta)
  }
  if (length(fixed) < 1) {
    m1 = VARorder(da, q + 12, output = FALSE)  # VARorfer  VAR模型的阶选择函数
    porder = m1$aicor                          # m1$aicor得到的是通过aic这个
    if (porder < 1) 
      porder = 1
    m2 = VAR(da, porder, output = FALSE)
    y = da[(porder + 1):nT, ]
    x = m2$residuals                           #residuals:残差
    m3 = THini(y, x, q, include.mean)          #函数VAR()中有:resi = y[, i] - xm %*% betai // 在这里残差就是正确值与预测值的差
    # 残差在数理统计中是指实际观察值与估计值(拟合值)之间的差。“残差”蕴含了有关模型基本假设的重要信息。如果回归模型正确的话, 我们可以将残差看作误差的观测值。
    beta = m3$estimates
    sebeta = m3$se
    nr = dim(beta)[1]
    if (prelim) {
      fixed = matrix(0, nr, k)
      for (j in 1:k) {
        tt = beta[, j]/sebeta[, j]
        idx = c(1:nr)[abs(tt) >= thres]
        fixed[idx, j] = 1
      }
    }
    if (length(fixed) < 1) {
      fixed = matrix(1, nr, k)
    }
  }
  else {
    nr = dim(beta)[1]
  }
  par = NULL
  separ = NULL
  fix1 = fixed
  VMAcnt = 0
  ist = 0
  if (include.mean) {
    jdx = c(1:k)[fix1[1, ] == 1]   #fixed中第一行中值是1的列数
    VMAcnt = length(jdx)
    if (VMAcnt > 0) {
      par = beta[1, jdx]
      separ = sebeta[1, jdx]
    }
    TH = -beta[2:(kq + 1), ]
    seTH = sebeta[2:(kq + 1), ]
    ist = 1
  }
  else {
    TH = -beta
    seTH = sebeta
  }
  for (j in 1:k) {
    idx = c(1:(nr - ist))[fix1[(ist + 1):nr, j] == 1]
    if (length(idx) > 0) {
      par = c(par, TH[idx, j])                  # 当fix为默认情况的时候,par=TH separ=seTH
      separ = c(separ, seTH[idx, j])
    }
  }
  ParMA <- par                                  # 当fix为默认情况的时候,ParMA=par=TH 
  LLKvma <- function(par, zt = zt, q = q, fixed = fix1, include.mean = include.mean) {
    k = ncol(zt)
    nT = nrow(zt)
    mu = rep(0, k)
    icnt = 0
    VMAcnt <- 0
    fix <- fixed
    iist = 0
    if (include.mean) {
      iist = 1
      jdx = c(1:k)[fix[1, ] == 1]
      icnt = length(jdx)
      VMAcnt <- icnt
      if (icnt > 0) 
        mu[jdx] = par[1:icnt]
    }
    for (j in 1:k) {
      zt[, j] = zt[, j] - mu[j]
    }
    kq = k * q
    Theta = matrix(0, kq, k)
    for (j in 1:k) {
      idx = c(1:kq)[fix[(iist + 1):(iist + kq), j] == 
                      1]
      jcnt = length(idx)
      if (jcnt > 0) {
        Theta[idx, j] = par[(icnt + 1):(icnt + jcnt)]
        icnt = icnt + jcnt
      }
    }
    TH = t(Theta)
    if (q > 1) {
      tmp = cbind(diag(rep(1, (q - 1) * k)), matrix(0, 
                                                    (q - 1) * k, k))
      TH = rbind(TH, tmp)
    }
    mm = eigen(TH)
    V1 = mm$values
    P1 = mm$vectors
    v1 = Mod(V1)
    ich = 0
    for (i in 1:kq) {
      if (v1[i] > 1) 
        V1[i] = 1/V1[i]
      ich = 1
    }
    if (ich > 0) {
      P1i = solve(P1)
      GG = diag(V1)
      TH = Re(P1 %*% GG %*% P1i)
      Theta = t(TH[1:k, ])
      ist = 0
      if (VMAcnt > 0) 
        ist = 1
      for (j in 1:k) {
        idx = c(1:kq)[fix[(ist + 1):(ist + kq), j] == 
                        1]
        jcnt = length(idx)
        if (jcnt > 0) {
          par[(icnt + 1):(icnt + jcnt)] = TH[j, idx]
          icnt = icnt + jcnt
        }
      }
    }
    at = mFilter(zt, t(Theta))
    sig = t(at) %*% at/nT
    ll = dmvnorm(at, rep(0, k), sig)
    LLKvma = -sum(log(ll))
    LLKvma
  }
  cat("Number of parameter s: ", length(par), "\n")
  cat("initial estimates: " , round(par, 4), "\n")
  lowerBounds = par
  upperBounds = par
  npar = length(par)
  mult = 2
  if ((npar > 10) || (q > 2)) 
    mult = 1.2
  for (j in 1:npar) {
    lowerBounds[j] = par[j] - mult * separ[j]
    upperBounds[j] = par[j] + mult * separ[j]
  }
  cat("Par. Lower-bounds: ", round(lowerBounds, 4), "\n")
  cat("Par. Upper-bounds: ", round(upperBounds, 4), "\n")
  if (details) {
    fit = nlminb(start = ParMA, objective = LLKvma, zt = da, 
                 fixed = fixed, include.mean = include.mean, q = q, 
                 lower = lowerBounds, upper = upperBounds, control = list(trace = 3))
  }
  else {
    fit = nlminb(start = ParMA, objective = LLKvma, zt = da, 
                 fixed = fixed, include.mean = include.mean, q = q, 
                 lower = lowerBounds, upper = upperBounds)
  }
  epsilon = 1e-04 * fit$par
  npar = length(par)
  Hessian = matrix(0, ncol = npar, nrow = npar)
  for (i in 1:npar) {
    for (j in 1:npar) {
      x1 = x2 = x3 = x4 = fit$par
      x1[i] = x1[i] + epsilon[i]
      x1[j] = x1[j] + epsilon[j]
      x2[i] = x2[i] + epsilon[i]
      x2[j] = x2[j] - epsilon[j]
      x3[i] = x3[i] - epsilon[i]
      x3[j] = x3[j] + epsilon[j]
      x4[i] = x4[i] - epsilon[i]
      x4[j] = x4[j] - epsilon[j]
      Hessian[i, j] = (LLKvma(x1, zt = da, q = q, fixed = fixed, 
                              include.mean = include.mean) - LLKvma(x2, zt = da, 
                                                                    q = q, fixed = fixed, include.mean = include.mean) - 
                         LLKvma(x3, zt = da, q = q, fixed = fixed, include.mean = include.mean) + 
                         LLKvma(x4, zt = da, q = q, fixed = fixed, include.mean = include.mean))/(4 * 
                                                                                                    epsilon[i] * epsilon[j])
    }
  }
  est = fit$par
  cat("Final   Estimates: ", est, "\n")
  se.coef = sqrt(diag(solve(Hessian)))
  tval = fit$par/se.coef
  matcoef = cbind(fit$par, se.coef, tval, 2 * (1 - pnorm(abs(tval))))
  dimnames(matcoef) = list(names(tval), c(" Estimate", " Std. Error", 
                                          " t value", "Pr(>|t|)"))
  cat("\nCoefficient(s):\n")
  printCoefmat(matcoef, digits = 4, signif.stars = TRUE)
  cat("---", "\n")
  cat("Estimates in matrix form:", "\n")
  icnt = 0
  ist = 0
  cnt = NULL
  if (include.mean) {
    ist = 1
    cnt = rep(0, k)
    secnt = rep(1, k)
    jdx = c(1:k)[fix1[1, ] == 1]
    icnt = length(jdx)
    if (icnt > 0) {
      cnt[jdx] = est[1:icnt]
      secnt[jdx] = se.coef[1:icnt]
      cat("Constant term: ", "\n")
      cat("Estimates: ", cnt, "\n")
    }
  }
  cat("MA coefficient matrix", "\n")
  TH = matrix(0, kq, k)
  seTH = matrix(1, kq, k)
  for (j in 1:k) {
    idx = c(1:kq)[fix1[(ist + 1):nr, j] == 1]
    jcnt = length(idx)
    if (jcnt > 0) {
      TH[idx, j] = est[(icnt + 1):(icnt + jcnt)]
      seTH[idx, j] = se.coef[(icnt + 1):(icnt + jcnt)]
      icnt = icnt + jcnt
    }
  }
  icnt = 0
  for (i in 1:q) {
    cat("MA(", i, ")-matrix", "\n")
    theta = t(TH[(icnt + 1):(icnt + k), ])
    print(theta, digits = 3)
    icnt = icnt + k
  }
  zt = da
  if (include.mean) {
    for (i in 1:k) {
      zt[, i] = zt[, i] - cnt[i]
    }
  }
  at = mFilter(zt, t(TH))
  sig = t(at) %*% at/nT
  cat(" ", "\n")
  cat("Residuals cov-matrix:", "\n")
  print(sig)
  dd = det(sig)
  d1 = log(dd)
  aic = d1 + 2 * npar/nT
  bic = d1 + log(nT) * npar/nT
  cat("----", "\n")
  cat("aic= ", aic, "\n")
  cat("bic= ", bic, "\n")
  Theta = t(TH)
  if (include.mean) {
    TH = rbind(cnt, TH)
    seTH = rbind(secnt, seTH)
  }
  VMA <- list(data = da, MAorder = q, cnst = include.mean, 
              coef = TH, secoef = seTH, residuals = at, Sigma = sig, 
              Theta = Theta, mu = cnt, aic = aic, bic = bic)
}

This snippet took 0.03 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).