# rates2.r # library(quadprog) # You may need to install the package # library(tseries) # You may need to install the package # rates <- read.table("rates.txt") rates <- as.matrix(rates) n=nrow(rates) year=as.matrix(rates[1:n,1]) month=as.matrix(rates[1:n,2]) t=year+(month-1)/12 r=as.matrix(rates[1:n,4]) dr=r[2:n]-r[1:(n-1)] n=n-1 # Create Data Matrix Using 24 Initial Conditions @ kk=24 # Number of initial conditions # nn=n-kk # Number of data points less number of initial conditions equals number of observations # y=as.matrix(dr[(1+kk):n]) # dependent variable # x=matrix(1,nn+1,1) # Regressors, first column (ones), one more observation than dependent variable (for forecast) # for (j in 1:kk) { x=cbind(x,dr[(1+kk-j):(n-j+1)]) } # X matrix columns, lags of y # # Model Combination # kn=kk+1 # Number of models = AR(0) through AR(kk) # yf=matrix(0,kn,1) # vector of forecasts (empty for now) # ee=matrix(0,nn,kn) # matrix of prediction errors (empty for now) # for (k in 1:kn) { xk=x[1:nn,1:k] xf=x[nn+1,1:k] xxi=solve(t(xk)%*%xk) beta=xxi%*%(t(xk)%*%y) e=y-xk%*%beta h=rowSums((xk%*%xxi)*xk) eh=e/(1-h) yf[k]=xf%*%beta ee[,k]=eh } Dmat=(t(ee)%*%ee)/nn dvec=matrix(0,kn,1) Amat=t(rbind(matrix(1,1,kn),diag(kn))) bvec=rbind(1,matrix(0,kn,1)) QP <- solve.QP(Dmat,dvec,Amat,bvec,bvec) w <- QP\$solution w <- as.matrix(w) e=ee%*%w cv=t(w)%*%Dmat%*%w yff=t(yf)%*%w print("Models, Weights") print(cbind(seq(0,kk),w)) print("CV, Combination Forecast, Level Forecast") print(cbind(cv,yff,r[n+1]+yff)) x.arch <- garch(e,order=c(1,1)) print(summary(x.arch)) archc=coef(x.arch) sd=predict(x.arch) like=logLik(x.arch) var <- as.matrix(sd[,1]^2) varf=archc[1]+archc[2]*(e[nn]^2)+archc[3]*var[nn,1] print("Forecast Variance, Standard Deviation") print(cbind(varf,sqrt(varf))) print("Unconditional Variance, Standard Deviation") sig=(t(e)%*%e)/nn print(cbind(sig,sqrt(sig))) print("Log Likelihood") print(like) t1 <- as.matrix(t[(kk+1):n]) e2=e^2 plot(t1,e,main="Leave-One-Out Prediction Residuals",type="l",xlab="",ylab="") windows() plot(t1,e2,main="Squared Prediction Residuals",type="l",xlab="",ylab="") windows() plot(t1,var,main="Estimated Variance",type="l",xlab="",ylab="")