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.