Dynamic Time Series Forecasting and Application

Question

I am trying to get a sliding dynamic time series prediction in R (and then work out the square prediction errors). I have a lot of this code based on https://stackoverflow.com/a/1386948/ ... , but I'm very new to R, so I'm struggling a bit. Any help is appreciated.

require(zoo) require(dynlm) set.seed(12345) #create variables x<-rnorm(mean=3,sd=2,100) y<-rep(NA,100) y[1]<-x[1] for(i in 2:100) y[i]=1+x[i-1]+0.5*y[i-1]+rnorm(1,0,0.5) int<-1:100 dummydata<-data.frame(int=int,x=x,y=y) zoodata<-as.zoo(dummydata) prediction<-function(series) { mod<-dynlm(formula = y ~ L(y) + L(x), data = series) #get model nextOb<-nrow(series)+1 #make forecast predicted<-coef(mod)[1]+coef(mod)[2]*zoodata$y[nextOb-1]+coef(mod)[3]*zoodata$x[nextOb-1] #strip timeseries information attributes(predicted)<-NULL return(predicted) } rolling<-rollapply(zoodata,width=40,FUN=prediction,by.column=FALSE) 

This returns:

 20 21 ..... 80 10.18676 10.18676 10.18676 

I have two problems that I did not expect:

• Works from 20-> 80, not 40-> 100, as you would expect (since the width is 40)
• The forecasts that he issues are constant: 10.18676

What am I doing wrong? And is there an easier way to make a prediction than to write all this? Thanks!

Answer 1

The main problem with your function is the data argument for dynlm . If you look in ?dynlm , you will see that the data argument must be a data.frame or zoo object. Unfortunately, I just found out that rollapply splits your zoo objects into array objects. This means that dynlm , noting that your data argument does not match the correct form, looked for x and y in your global environment, which, of course, were defined at the top of your code. The solution is to convert the series to a zoo object. There were a few more problems in your code, I posted the corrected version here:

 prediction<-function(series) { mod <- dynlm(formula = y ~ L(y) + L(x), data = as.zoo(series)) # get model # nextOb <- nrow(series)+1 # This will always be 21. I think you mean: nextOb <- max(series[,'int'])+1 # To get the first row that follows the window if (nextOb<=nrow(zoodata)) { # You won't predict the last one # make forecast # predicted<-coef(mod)[1]+coef(mod)[2]*zoodata$y[nextOb-1]+coef(mod)[3]*zoodata$x[nextOb-1] # That would work, but there is a very nice function called predict predicted=predict(mod,newdata=data.frame(x=zoodata[nextOb,'x'],y=zoodata[nextOb,'y'])) # I'm not sure why you used nextOb-1 attributes(predicted)<-NULL # I added the square error as well as the prediction. c(predicted=predicted,square.res=(predicted-zoodata[nextOb,'y'])^2) } } rollapply(zoodata,width=20,FUN=prediction,by.column=F,align='right') 

Your second question about numbering your results can be controlled by the align rollapply argument. left will give you 1..60 , center (by default) will give you 20..80 and right will get you 40..100 .