Recursive forecasting in R

Question

setwd("C:/Users/user/Desktop")
library(forecast)
library(vars)
library(tseries)
cay <- read.csv("C:/Users/User/Desktop/cay.csv")
sp500 <- read.csv("C:/Users/User/Desktop/sp500.csv")
cay_indicator <- cay$indicator
cay <- cay$cay
sp500c <- sp500[,2]
sp500c <- data.matrix(sp500c)
sp500c.ts=ts(data = sp500c, start = c(1950,2), end = c(2015,4), frequency = 4)
plot(sp500c.ts)
sp500=ts(data = sp500c, start = c(1953, 1), end = c(2014, 3), frequency = 4)
cay=ts(data = cay, start = c(1953, 2), end = c(2014, 3), frequency = 4)
pp.test(sp500)
adf.test(sp500)
returnsp500=diff(log(sp500))
pp.test(returnsp500)
adf.test(returnsp500)
var1=ts(cbind(returnsp500, cay), start = c(1953, 2), end = c(2014, 3), frequency = 4)
var1     
var2=ts(cbind(returnsp500, cay_indicator), start = c(1953, 2), end = c(2014, 3), frequency = 4)
var2
var3=ts(cbind(returnsp500, cay, cay_indicator), start = c(1953, 2), end = c(2014, 3), frequency = 4)
var3
layout(1)
plot(var3,type="l",col="blue",main="Stock Price, Consumption Wealth Ratio, Negativity Indicator")
acf(var1, lag.max=24)
acf(var2, lag.max=24)
acf(var3,lag.max=24)
info.crit1=VARselect(var1,lag.max=5,type="const")
info.crit2=VARselect(var2,lag.max=5,type="const")
info.crit3=VARselect(var3,lag.max=5,type="const")
info.crit1
info.crit2
info.crit3
model1=VAR(var1,p=1,type="const")
model1
summary(model1)
model2=VAR(var2,p=1,type="const")
model2
summary(model2)
model3=VAR(var3,p=1,type="const")
model3
summary(model3)
causality(model1,cause="cay")$Granger
causality(model2,cause="cay_indicator")$Granger

n.end=127
t=length(returnsp500)
n=t-n.end-3

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=1,se.fit=FALSE)
  pred1[i,1]=for_var$fcst$returnsp500[,1]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=1,se.fit=FALSE)
  pred2[i,1]=for_var$fcst$returnsp500[,1]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=1,se.fit=FALSE)
  pred3[i,1]=for_var$fcst$returnsp500[,1]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1]=predict(model.ar,n.ahead=1,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,1),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,1),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,1),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,1),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+1):(t-3)],start=c(1985,1),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,1],col="green")
lines(pred2.ts[,1],col="red")
lines(pred3.ts[,1],col="black")
lines(pred4.ts[,1],col="yellow")

pred1.ts

e1_var=returnsp500[(n.end+1):(t-3)]-pred1.ts[,1] #1-step ahead bivariate cay & sp500
e2_var=returnsp500[(n.end+1):(t-3)]-pred2.ts[,1] #1-step ahead bivariate cay_indicator & sp500
e3_var=returnsp500[(n.end+1):(t-3)]-pred3.ts[,1] #1-step ahead bivariate cay and cay_indicator & sp500
e4_var=returnsp500[(n.end+1):(t-3)]-pred4.ts[,1] #1-step ahead univariate model just a lagged sp500 

n.end=127
t=length(returnsp500) 
n=t-n.end-3

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=2,se.fit=FALSE)
  pred1[i,1:2]=for_var$fcst$returnsp500[1:2]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=2,se.fit=FALSE)
  pred2[i,1:2]=for_var$fcst$returnsp500[1:2]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=2,se.fit=FALSE)
  pred3[i,1:2]=for_var$fcst$returnsp500[1:2]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1:2]=predict(model.ar,n.ahead=2,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,2),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,2),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,2),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,2),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+2):(t-2)],start=c(1985,2),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,2],col="green")
lines(pred2.ts[,2],col="red")
lines(pred3.ts[,2],col="black")
lines(pred4.ts[,2],col="yellow")


e1_var2=returnsp500[(n.end+2):(t-2)]-pred1.ts[,2] #2-step ahead bivariate cay & sp500
e2_var2=returnsp500[(n.end+2):(t-2)]-pred2.ts[,2] #2-step ahead bivariate cay_indicator & sp500
e3_var2=returnsp500[(n.end+2):(t-2)]-pred3.ts[,2] #2-step ahead bivariate cay and cay_indicator & sp500
e4_var2=returnsp500[(n.end+2):(t-2)]-pred4.ts[,2] #2-step ahead univariate model just a lagged sp500 


n.end=127
t=length(returnsp500) 
n=t-n.end-3

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=3,se.fit=FALSE)
  pred1[i,1:3]=for_var$fcst$returnsp500[1:3]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=3,se.fit=FALSE)
  pred2[i,1:3]=for_var$fcst$returnsp500[1:3]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=3,se.fit=FALSE)
  pred3[i,1:3]=for_var$fcst$returnsp500[1:3]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1:3]=predict(model.ar,n.ahead=3,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,3),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,3),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,3),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,3),frequency=4)
logsp500returns.ts=ts(data=returnsp500,start=c(1985,3), end = c(2014,2),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,3],col="green")
lines(pred2.ts[,3],col="red")
lines(pred3.ts[,3],col="black")
lines(pred4.ts[,3],col="yellow")



e1_var3=returnsp500[(n.end+2):(t-2)]-pred1.ts[,3] #3-step ahead bivariate cay & sp500
e2_var3=returnsp500[(n.end+2):(t-2)]-pred2.ts[,3] #3-step ahead bivariate cay_indicator & sp500
e3_var3=returnsp500[(n.end+2):(t-2)]-pred3.ts[,3] #3-step ahead bivariate cay and cay_indicator & sp500
e4_var3=returnsp500[(n.end+2):(t-2)]-pred4.ts[,3] #3-step ahead univariate model just a lagged sp500 


n.end=127
t=length(returnsp500) 
n=t-n.end-3

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=4,se.fit=FALSE)
  pred1[i,1:4]=for_var$fcst$returnsp500[1:4]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=4,se.fit=FALSE)
  pred2[i,1:4]=for_var$fcst$returnsp500[1:4]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=4,se.fit=FALSE)
  pred3[i,1:4]=for_var$fcst$returnsp500[1:4]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1:4]=predict(model.ar,n.ahead=4,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,4),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,4),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,4),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,4),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+4):t],start=c(1985,4),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,4],col="green")
lines(pred2.ts[,4],col="red")
lines(pred3.ts[,4],col="black")
lines(pred4.ts[,4],col="yellow")



e1_var4=returnsp500[(n.end+4):(t)]-pred1.ts[,4] #4-step ahead bivariate cay & sp500
e2_var4=returnsp500[(n.end+4):(t)]-pred2.ts[,4] #4-step ahead bivariate cay_indicator & sp500
e3_var4=returnsp500[(n.end+4):(t)]-pred3.ts[,4] #4-step ahead bivariate cay and cay_indicator & sp500
e4_var4=returnsp500[(n.end+4):(t)]-pred4.ts[,4] #4-step ahead univariate model just a lagged sp500 

rmse1_var1=sqrt(mean(e1_var^2)) #16. RMSE's in order, 1-step, 2-step, 3-step, 4-step
rmse2_var1=sqrt(mean(e2_var^2))
rmse3_var1=sqrt(mean(e3_var^2))
rmse4_var1=sqrt(mean(e4_var^2))

rmse1_var2=sqrt(mean(e1_var2^2)) 
rmse2_var2=sqrt(mean(e2_var2^2))
rmse3_var2=sqrt(mean(e3_var2^2))
rmse4_var2=sqrt(mean(e4_var2^2))

rmse1_var3=sqrt(mean(e1_var3^2)) 
rmse2_var3=sqrt(mean(e2_var3^2))
rmse3_var3=sqrt(mean(e3_var3^2))
rmse4_var3=sqrt(mean(e4_var3^2))

rmse1_var4=sqrt(mean(e1_var4^2)) 
rmse2_var4=sqrt(mean(e2_var4^2))
rmse3_var4=sqrt(mean(e3_var4^2))
rmse4_var4=sqrt(mean(e4_var4^2))

rmse1_var1
rmse2_var1                                    
rmse3_var1
rmse4_var1

rmse1_var2
rmse2_var2
rmse3_var2
rmse4_var2

rmse1_var3
rmse2_var3
rmse3_var3
rmse4_var3

rmse1_var4
rmse2_var4
rmse3_var4
rmse4_var4

c=rep(1,116)
cw=e4_var^2-e1_var^2+(e4_var-e1_var)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var2^2-e1_var2^2+(e4_var2-e1_var2)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var3^2-e1_var3^2+(e4_var3-e1_var3)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var4^2-e1_var4^2+(e4_var4-e1_var4)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var^2-e2_var^2+(e4_var-e2_var)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var2^2-e2_var2^2+(e4_var2-e2_var2)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var3^2-e2_var3^2+(e4_var3-e2_var3)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var4^2-e2_var4^2+(e4_var4-e2_var4)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

dmw=(e2_var^2-e1_var^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

dmw=(e2_var2^2-e1_var2^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

dmw=(e2_var3^2-e1_var3^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

dmw=(e2_var4^2-e1_var4^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

c=rep(1,116)
cw=e1_var^2-e3_var^2+(e1_var-e3_var)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e1_var2^2-e3_var2^2+(e1_var2-e3_var2)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e1_var3^2-e3_var3^2+(e1_var3-e3_var3)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e1_var4^2-e3_var4^2+(e1_var4-e3_var4)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

Here is my r code, the for loop is my only problem seemingly. I am using 3 var models, 2 bivariates and a trivariate. I need help getting the for loop to correctly give me predictions. I cannot be certain that it is doing this right now, in fact I have a hunch that it is not. If someone could please tell me in my for loop what is incorrect, it would be much appreciated, or at least how to fix my for loop.

You can also shorten your code quite a bit by writing functions for all the copy and paste that you did — Tung, Apr 12 '18 at 04:00
It's hard to tell, try giving a smaller example of where you think it is not right. Also we do not have your data. To debug reduce your loop to one model and incremente complexity as you check your code to make sure it is doing what it is supposed to. — R. Prost, Apr 12 '18 at 05:19
I changed my code, it's still quite long, so pardon me for that. But this might be it actually. Just for reference, R Prost you did a good job of giving some quick but useful help to help me focus on what I was doing incorrectly. Any help is still welcomed. — Nicolas Cameranesi, Apr 12 '18 at 06:11

score 0 · Answer 1 · answered Apr 12 '18 at 06:12

This is a possible ans.

setwd("C:/Users/user/Desktop")
library(forecast)
library(vars)
library(tseries)
cay <- read.csv("C:/Users/User/Desktop/cay.csv")
sp500 <- read.csv("C:/Users/User/Desktop/sp500.csv")
cay_indicator <- cay$indicator
cay <- cay$cay
sp500c <- sp500[,2]
sp500c <- data.matrix(sp500c)
sp500c.ts=ts(data = sp500c, start = c(1950,2), end = c(2015,4), frequency = 4)
plot(sp500c.ts)
sp500=ts(data = sp500c, start = c(1953, 1), end = c(2014, 3), frequency = 4)
cay=ts(data = cay, start = c(1953, 2), end = c(2014, 3), frequency = 4)
pp.test(sp500)
adf.test(sp500)
returnsp500=diff(log(sp500))
pp.test(returnsp500)
adf.test(returnsp500)
var1=ts(cbind(returnsp500, cay), start = c(1953, 2), end = c(2014, 3), frequency = 4)
var1     
var2=ts(cbind(returnsp500, cay_indicator), start = c(1953, 2), end = c(2014, 3), frequency = 4)
var2
var3=ts(cbind(returnsp500, cay, cay_indicator), start = c(1953, 2), end = c(2014, 3), frequency = 4)
var3
layout(1)
plot(var3,type="l",col="blue",main="Stock Price, Consumption Wealth Ratio, Negativity Indicator")
acf(var1, lag.max=24)
acf(var2, lag.max=24)
acf(var3,lag.max=24)
info.crit1=VARselect(var1,lag.max=5,type="const")
info.crit2=VARselect(var2,lag.max=5,type="const")
info.crit3=VARselect(var3,lag.max=5,type="const")
info.crit1
info.crit2
info.crit3
model1=VAR(var1,p=1,type="const")
model1
summary(model1)
model2=VAR(var2,p=1,type="const")
model2
summary(model2)
model3=VAR(var3,p=1,type="const")
model3
summary(model3)
causality(model1,cause="cay")$Granger
causality(model2,cause="cay_indicator")$Granger

n.end=127
t=length(returnsp500)
n=t-n.end

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=1,se.fit=FALSE)
  pred1[i,1]=for_var$fcst$returnsp500[,1]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=1,se.fit=FALSE)
  pred2[i,1]=for_var$fcst$returnsp500[,1]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=1,se.fit=FALSE)
  pred3[i,1]=for_var$fcst$returnsp500[,1]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1]=predict(model.ar,n.ahead=1,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,1),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,1),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,1),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,1),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+1):(t)],start=c(1985,1),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,1],col="green")
lines(pred2.ts[,1],col="red")
lines(pred3.ts[,1],col="black")
lines(pred4.ts[,1],col="yellow")

pred1.ts[,1]

e1_var=returnsp500[(n.end+1):(t)]-pred1.ts[,1] #1-step ahead bivariate cay & sp500
e2_var=returnsp500[(n.end+1):(t)]-pred2.ts[,1] #1-step ahead bivariate cay_indicator & sp500
e3_var=returnsp500[(n.end+1):(t)]-pred3.ts[,1] #1-step ahead bivariate cay and cay_indicator & sp500
e4_var=returnsp500[(n.end+1):(t)]-pred4.ts[,1] #1-step ahead univariate model just a lagged sp500 

n.end=127
t=length(returnsp500) 
n=t-n.end-1

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=2,se.fit=FALSE)
  pred1[i,1:2]=for_var$fcst$returnsp500[1:2]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=2,se.fit=FALSE)
  pred2[i,1:2]=for_var$fcst$returnsp500[1:2]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=2,se.fit=FALSE)
  pred3[i,1:2]=for_var$fcst$returnsp500[1:2]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1:2]=predict(model.ar,n.ahead=2,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,2),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,2),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,2),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,2),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+2):(t)],start=c(1985,2),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,2],col="green")
lines(pred2.ts[,2],col="red")
lines(pred3.ts[,2],col="black")
lines(pred4.ts[,2],col="yellow")

pred1.ts[,2]

e1_var2=returnsp500[(n.end+2):(t)]-pred1.ts[,2] #2-step ahead bivariate cay & sp500
e2_var2=returnsp500[(n.end+2):(t)]-pred2.ts[,2] #2-step ahead bivariate cay_indicator & sp500
e3_var2=returnsp500[(n.end+2):(t)]-pred3.ts[,2] #2-step ahead bivariate cay and cay_indicator & sp500
e4_var2=returnsp500[(n.end+2):(t)]-pred4.ts[,2] #2-step ahead univariate model just a lagged sp500 


n.end=127
t=length(returnsp500) 
n=t-n.end-2

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=3,se.fit=FALSE)
  pred1[i,1:3]=for_var$fcst$returnsp500[1:3]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=3,se.fit=FALSE)
  pred2[i,1:3]=for_var$fcst$returnsp500[1:3]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=3,se.fit=FALSE)
  pred3[i,1:3]=for_var$fcst$returnsp500[1:3]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1:3]=predict(model.ar,n.ahead=3,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,3),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,3),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,3),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,3),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+1):(t)],start=c(1985,3), frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,3],col="green")
lines(pred2.ts[,3],col="red")
lines(pred3.ts[,3],col="black")
lines(pred4.ts[,3],col="yellow")

e1_var3=returnsp500[(n.end+3):(t)]-pred1.ts[,3] #3-step ahead bivariate cay & sp500
e2_var3=returnsp500[(n.end+3):(t)]-pred2.ts[,3] #3-step ahead bivariate cay_indicator & sp500
e3_var3=returnsp500[(n.end+3):(t)]-pred3.ts[,3] #3-step ahead bivariate cay and cay_indicator & sp500
e4_var3=returnsp500[(n.end+3):(t)]-pred4.ts[,3] #3-step ahead univariate model just a lagged sp500 


n.end=127
t=length(returnsp500) 
n=t-n.end-3

pred1=matrix(rep(0,4*n),n,4)
pred2=matrix(rep(0,4*n),n,4)
pred3=matrix(rep(0,4*n),n,4)
pred4=matrix(rep(0,4*n),n,4)

for(i in 1:n){
  x_var=var1[1:n.end+i-1,]
  x_var2=var2[1:n.end+i-1,]
  x_var3=var3[1:n.end+i-1,]
  x_ar=returnsp500[1:n.end+i-1]

  model.var=VAR(x_var,p=1,type="const")
  for_var=predict(model.var,n.ahead=4,se.fit=FALSE)
  pred1[i,1:4]=for_var$fcst$returnsp500[1:4]

  model.var=VAR(x_var2,p=1,type="const")
  for_var=predict(model.var,n.ahead=4,se.fit=FALSE)
  pred2[i,1:4]=for_var$fcst$returnsp500[1:4]

  model.var=VAR(x_var3,p=1,type="const")
  for_var=predict(model.var,n.ahead=4,se.fit=FALSE)
  pred3[i,1:4]=for_var$fcst$returnsp500[1:4]

  model.ar=arima(x_ar,order=c(1,0,0),method="ML")
  pred4[i,1:4]=predict(model.ar,n.ahead=4,se.fit=FALSE)[1]
}

pred1.ts=ts(data=pred1,start=c(1985,4),frequency=4)
pred2.ts=ts(data=pred2,start=c(1985,4),frequency=4)
pred3.ts=ts(data=pred3,start=c(1985,4),frequency=4)
pred4.ts=ts(data=pred4,start=c(1985,4),frequency=4)
logsp500returns.ts=ts(data=returnsp500[(n.end+4):t],start=c(1985,4),frequency=4)
plot(logsp500returns.ts,col="blue",ylim=c(-.1,.1))
lines(pred1.ts[,4],col="green")
lines(pred2.ts[,4],col="red")
lines(pred3.ts[,4],col="black")
lines(pred4.ts[,4],col="yellow")


e1_var4=returnsp500[(n.end+4):(t)]-pred1.ts[,4] #4-step ahead bivariate cay & sp500
e2_var4=returnsp500[(n.end+4):(t)]-pred2.ts[,4] #4-step ahead bivariate cay_indicator & sp500
e3_var4=returnsp500[(n.end+4):(t)]-pred3.ts[,4] #4-step ahead bivariate cay and cay_indicator & sp500
e4_var4=returnsp500[(n.end+4):(t)]-pred4.ts[,4] #4-step ahead univariate model just a lagged sp500 

rmse1_var1=sqrt(mean(e1_var^2)) #16. RMSE's in order, 1-step, 2-step, 3-step, 4-step
rmse2_var1=sqrt(mean(e2_var^2))
rmse3_var1=sqrt(mean(e3_var^2))
rmse4_var1=sqrt(mean(e4_var^2))

rmse1_var2=sqrt(mean(e1_var2^2)) 
rmse2_var2=sqrt(mean(e2_var2^2))
rmse3_var2=sqrt(mean(e3_var2^2))
rmse4_var2=sqrt(mean(e4_var2^2))

rmse1_var3=sqrt(mean(e1_var3^2)) 
rmse2_var3=sqrt(mean(e2_var3^2))
rmse3_var3=sqrt(mean(e3_var3^2))
rmse4_var3=sqrt(mean(e4_var3^2))

rmse1_var4=sqrt(mean(e1_var4^2)) 
rmse2_var4=sqrt(mean(e2_var4^2))
rmse3_var4=sqrt(mean(e3_var4^2))
rmse4_var4=sqrt(mean(e4_var4^2))

rmse1_var1
rmse2_var1                                    
rmse3_var1
rmse4_var1

rmse1_var2
rmse2_var2
rmse3_var2
rmse4_var2

rmse1_var3
rmse2_var3
rmse3_var3
rmse4_var3

rmse1_var4
rmse2_var4
rmse3_var4
rmse4_var4

c=rep(1,116)
cw=e4_var^2-e1_var^2+(e4_var-e1_var)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var2^2-e1_var2^2+(e4_var2-e1_var2)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var3^2-e1_var3^2+(e4_var3-e1_var3)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var4^2-e1_var4^2+(e4_var4-e1_var4)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var^2-e2_var^2+(e4_var-e2_var)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var2^2-e2_var2^2+(e4_var2-e2_var2)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var3^2-e2_var3^2+(e4_var3-e2_var3)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e4_var4^2-e2_var4^2+(e4_var4-e2_var4)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

dmw=(e2_var^2-e1_var^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

dmw=(e2_var2^2-e1_var2^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

dmw=(e2_var3^2-e1_var3^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

dmw=(e2_var4^2-e1_var4^2)
c=rep(1,116)
reg=lm(dmw~c-1)
avar=NeweyWest(reg,lag=NULL,prewhite=FALSE)
dmw.q=reg$coef/sqrt(avar)
dmw.q
pnorm(dmw.q)
p.value=1-pnorm(dmw.q)
pstr.value=p.value/2 
pstr.value

c=rep(1,116)
cw=e1_var^2-e3_var^2+(e1_var-e3_var)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e1_var2^2-e3_var2^2+(e1_var2-e3_var2)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e1_var3^2-e3_var3^2+(e1_var3-e3_var3)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

c=rep(1,116)
cw=e1_var4^2-e3_var4^2+(e1_var4-e3_var4)^2
reg.cw=lm(cw~c-1)
avar.cw=NeweyWest(reg.cw,lag=NULL,prewhite=FALSE)
cw.test=reg.cw$coef/sqrt(avar.cw)
cw.test
pnorm(cw.test)
(1-pnorm(cw.test))/2

Recursive forecasting in R

1 Answers1