jpm miao
2012-Aug-27 08:02 UTC
[R] How can I find the principal components and run regression/forecasting using dynlm
Hello,
I would like to write a program that compute the principal components of
a set of data and then
1. Run the dependent variable against the principal components (lagged
value)
2. Do prediction
, following Stock and Watson (1999) "Forecasting Inflation". All
data
are time series.
Now I can run the program using single factor (first principal
component), but I don't know how to write a program for muli-factor model
(first several principal component).
The core subroutine called is dynlm, which runs OLS for time series. I
think the primary problem is that (1) I don't how to run dynlm with matrix
rather than vector as explanatory variables (2)I don't know how to do a
forecast with the estimation results of dynlm properly. In lm model,
function "predict.lm" can do it.
For the case of first principal component (In order to accentuate the
main problem, I have simplify the codes):
Mpca<-prcomp(M1, center=TRUE, scale =TRUE) # M1 is the data matrix of
explanatory variables
Mpca1st<-Mpca$x[,1] # first principle component
X<-as.matrix(Mpca1st)
model<-dynlm(as.ts(y[(h+1):t]) ~ L(as.ts(X[1:(t-h)]), 0:i) +
L(as.ts(z[1:(t-h)]),0:j)) # y, X, z are a zoo objects defined earlier; i
and j, t, h are given earlier; h is the forecasting horizon; since dynlm
can't work for zoo object, the variables need to be transformed to ts
objects.
Xlast<-as.matrix(X[(t-h):(t-h-i+1)])
zlast<-t(as.matrix(z[(t-h):(t-h-j+1)])) # Make it a column matrix
beta<-as.matrix(model$coefficient)
pi1fore[t]<-pi1[t-h]+c(1, Xlast, zIlast)%*%beta # pi1fore is defined
earlier . I just use the inner product as a vehicle for forecasting. I
don't know if there's a way that the codes can be written neatly?
return(pi1fore)
Then I attempt to program similarly for the multi-factor case:
Mpca<-prcomp(M1, center=TRUE, scale =TRUE) # M1 is the data matrix of
explanatory variables
Mpca1f<-Mpca$x[,(1:f)] # first f principle components
X<-as.matrix(Mpca1f)
model<-dynlm(as.ts(y[(h+1):t]) ~ L(as.ts(X[1:(t-h),]), 0:i) +
L(as.ts(z[1:(t-h)]),0:j)) # Since X has more than one column, I add a
comma after 1:(t-h). I don't know if I am right here?
Xlast<-as.matrix(X[(t-h):(t-h-i+1), ])
zlast<-t(as.matrix(z[(t-h):(t-h-j+1)])) # Make it a column matrix
beta<-as.matrix(model$coefficient)
pi1fore[t]<-pi1[t-h]+c(1, Xlast, zIlast)%*%beta # It does not seem to
work here! How can I modify the program...?.
return(pi1fore)
[[alternative HTML version deleted]]
