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]]