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

这是我的 r 代码，for 循环似乎是我唯一的问题。我正在使用 3 个 var 模型、2 个二元变量和一个三元变量。我需要帮助让 for 循环正确地给我预测。我不能确定它现在正在这样做，事实上我有一种预感，它不是。如果有人可以在我的 for 循环中告诉我什么是不正确的，将不胜感激，或者至少如何修复我的 for 循环。

Answer 1

这是一个可能的答案。

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