8
$\begingroup$

I would like to predict using a linear model in R. The model that I have is of the form lm(y~ lag(x)). It seems like I should be able to predict using the predict function, but this does not look ahead into the future.

Here is the lag function that I have been using from within R.

lag1 = function (x) c(NA, x[1:(length(x)-1)])

This lag function adds an NA at the beginning of the data, shifts everything down, and cuts off the last observation.

I am mostly interested in the predictions made on the last line of the data. This line of data seems to be ignored by the fitted and predict functions.

The best workaround I have found was to multiply each of the model coefficients by each of my dependent variables, and add them up for each prediction. None of the dependent variables are lagged at this point, when I set up this equation. The problem with this is that each time I change the model formula, the names of the variables have to be changed to match the new formula. It seems like there should be a more intuitive solution to this.

$\endgroup$
3
  • $\begingroup$ Not sure I understand. Did you give new x values to the predict() function as the second argument? $\endgroup$ Commented Apr 4, 2014 at 6:37
  • 1
    $\begingroup$ Simply feed the future regressor values to the newdata parameter of predict.lm(). Make sure they are in a data.frame with the same names as used in the formula in lm(). $\endgroup$ Commented Apr 4, 2014 at 6:51
  • $\begingroup$ @Stephan Kollasa, For time series data your solution unfortunately does not work. $\endgroup$
    – mpiktas
    Commented Apr 4, 2014 at 8:20

1 Answer 1

19
$\begingroup$

Estimating with lags and using model for predicting is a sore point in base R. Here is the example:

set.seed(1)
y<-ts(rnorm(10))
x<-ts(rnorm(10))
lm(y~x)

Call:
lm(formula = y ~ x)

Coefficients:
(Intercept)            x  
     0.2006      -0.2749  

>    lm(y~lag(x))

Call:
lm(formula = y ~ lag(x))

Coefficients:
(Intercept)       lag(x)  
     0.2006      -0.2749  

Notice that in both cases the result is the same. This happens because of the somewhat peculiar behaviour of the default lag function. Compare x and lag(x)

 x
Time Series:
Start = 1 
End = 10 
Frequency = 1 
 [1]  1.51178117  0.38984324 -0.62124058 -2.21469989  1.12493092 -0.04493361 -0.01619026  0.94383621  0.82122120
[10]  0.59390132
> lag(x)
Time Series:
Start = 0 
End = 9 
Frequency = 1 
 [1]  1.51178117  0.38984324 -0.62124058 -2.21469989  1.12493092 -0.04493361 -0.01619026  0.94383621  0.82122120
[10]  0.59390132

As you can see the data is the same, only the attributes, in this case the time, are different. Hence the lm sees the same data, since it ignores the attributes. There are several ways of working around this behaviour. Here are few.

First you can convert the data to time series format for which the lag behaviour is "standard". One of such formats is xts from package xts:

yy<-xts(y,as.Date(1:10))
xx<-xts(x,as.Date(1:10))
lm(yy~xx)

Call:
lm(formula = yy ~ lag(xx))

Coefficients:
(Intercept)      lag(xx)  
     0.2754      -0.2798 

See now the coefficient is different, since lag now correctly shifts the data:

lag(xx)
                  [,1]
1970-01-02          NA
1970-01-03  1.51178117
1970-01-04  0.38984324
1970-01-05 -0.62124058
1970-01-06 -2.21469989
1970-01-07  1.12493092
1970-01-08 -0.04493361
1970-01-09 -0.01619026
1970-01-10  0.94383621
1970-01-11  0.82122120

Another way is to estimate the regression using the function dynlm from package dynlm:

dynlm(y~L(x))

Time series regression with "ts" data:
Start = 2, End = 10

Call:
dynlm(formula = y ~ L(x))

Coefficients:
(Intercept)         L(x)  
     0.2754      -0.2798  

This covers the estimation. Now predicting is trickier. To get the fitted values you can simply use predict function:

predict(dynlm(y ~ L(x)))
          2           3           4           5           6           7           8           9          10 
-0.14757748  0.16632219  0.44920672  0.89503027 -0.03934318  0.28796556  0.27992365  0.01132412  0.04562977 
predict(lm(yy ~ lag(xx)))
          2           3           4           5           6           7           8           9          10 
-0.14757748  0.16632219  0.44920672  0.89503027 -0.03934318  0.28796556  0.27992365  0.01132412  0.04562977 

Predicting ahead in a future though presents a problem. The default behaviour of predict function is to expect newdata argument. But for one step ahead forecast in this case, no new data is required. So standard predict function will not work in this case. I would love to see a general solution for this problem, but in my knowledge different packages provide different ways of getting such forecasts and if you do not want to write the function yourself you need to cast your model in a form which is required by a specific forecasting function from a specific package. And you need to know that package pretty well.

One of such packages is midasr (of which I am developer). One step ahead forecast for such model would be implemented in the following way:

midas_u(y~mls(x,1,1))

Call:
lm(formula = y ~ mls(x, 1, 1), data = ee)

Coefficients:
 (Intercept)  mls(x, 1, 1)  
      0.2754       -0.2798  

forecast(midas_u(y~mls(x,1,1)),newdata=list(x=NA))
[1] 0.1092301

The package midasr works with mixed frequency data. The function mls has 3 arguments, data, lag numbers and frequency ratio. In this case the frequency is the same so the third argument is 1. In that case the function works exactly as the function lag for xts objects. For forecasting it is necessary to supply new data. As we do not know it we supply the NA which works in this case, since one-step ahead forecasts in this case needs only the data which we already know.

The last value of x is 0.59390132, so you can check that the result is the correct one directly

> 0.59390132*(-0.2798)  +0.2754
[1] 0.1092264

The answer is correct to 4 decimal places, since I used the coefficients with 4 digit accuracy.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.