/
forecast.R
105 lines (73 loc) · 2.42 KB
/
forecast.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
#' @export
#' @title forecasting using an ltvar model
#' @param obj a fitted ltvar model
#' @param n.ahead forecast horizon
#' @param ... currently not used
forecast <- function(obj,n.ahead,...) UseMethod("forecast")
#' @export
#' @title forecasting using an ltvar model
#' @param obj a fitted ltvar model
#' @param n.ahead forecast horizon
#' @param ... currently not used
#' @importFrom stats rnorm
forecast.ltvar <- function(obj,n.ahead,...){
# Preliminaries
y = obj$data_info$data
p = obj$model_info$p
intercept = obj$model_info$intercept
nreps <- dim(obj$mcmc_draws$mb)[1]
nTotal <- dim(obj$mcmc_draws$mb)[2]
nb <- dim(obj$mcmc_draws$mb)[3]
na <- dim(obj$mcmc_draws$ma)[3]
nk <- dim(obj$mcmc_draws$mh)[3]
forecast_storage <- array(0,dim=c(n.ahead,nk,nreps))
for(ii in 1:nreps){
mb <- obj$mcmc_draws$mb[ii,nTotal,]
vmub <- obj$mcmc_draws$vmub[ii,,]
mSigb <- obj$mcmc_draws$mSigb[ii,,]
mPhib <- obj$mcmc_draws$mPhib[ii,,]
vdb <- obj$mcmc_draws$vdb[ii,,]
ma <- obj$mcmc_draws$ma[ii,nTotal,]
vmua <- obj$mcmc_draws$vmua[ii,,]
mSiga <- obj$mcmc_draws$mSiga[ii,,]
mPhia <- obj$mcmc_draws$mPhia[ii,,]
vda <- obj$mcmc_draws$vda[ii,,]
mPhih <- obj$mcmc_draws$mPhih[ii,,]
mSigh <- obj$mcmc_draws$mSigh[ii,,]
vbf <- obj$mcmc_draws$mb[ii, nTotal - 1,]
vaf <- obj$mcmc_draws$ma[ii, nTotal - 1,]
vhf <- obj$mcmc_draws$mh[ii, nTotal - 1,]
vgam <- obj$mcmc_draws$vgam[ii,,]
myf <- y[nTotal:(nTotal - p + 1),]
for(jj in 1:(n.ahead)){
# Calculate time-varying parameters
vbf <- vmub + mPhib %*% (vbf - vmub) + t(chol(mSigb)) %*% stats::rnorm(nb)
vbfd <- vbf
vbfd[abs(vbf) < vdb] <- 0
vaf <- vmua + mPhia %*% (vaf - vmua) + t(chol(mSiga)) %*% stats::rnorm(na)
vafd <- vaf
vafd[abs(vaf) < vda] <- 0
vhf <- mPhih %*% vhf + t(chol(mSigh)) %*% rnorm(nk)
mOmsq <- solve(fAt(vafd,nk)) %*% diag(exp(as.vector(vhf)/2)) * sqrt(vgam)
vb <- matrix(vbfd,ncol = nk)
if(intercept){
tmp <-c(1,matrix(t(myf),nrow = 1))
}
else{
tmp <- matrix(t(myf),nrow = 1)
}
vyf <- t(tmp %*% vb) + t(chol(mOmsq)) %*% rnorm(nk)
# Geet new variables
if(p > 1){
myf[2:p,] <- myf[1:(p-1),]
myf[1,] <- vyf
}
else{
myf <- vyf
}
# Store forecasts
forecast_storage[jj,,ii] <- vyf
}
}
return(forecast_storage)
}