# rates1.r # 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]) plot(t,r,main="10-Year Treasury Bill Rate",type="l",xlab="",ylab="") # First Difference # dr=r[2:n]-r[1:(n-1)] windows() plot(t[2:n],dr,main="Change in 10-Year Rate",type="l",xlab="",ylab="") 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 # # Estimate AR(12) # k=13 # Number of lags + 1 for intercept # xk=x[1:nn,1:k] # X matrix, for estimation # xf=x[nn+1,1:k] # X vector, for forecast # xxi=solve(t(xk)%*%xk) # X'X # beta=xxi%*%(t(xk)%*%y) # OLS # e=y-xk%*%beta # OLS residual # h=rowSums((xk%*%xxi)*xk) # leverage values # eh=e/(1-h) # LOO prediction residuals # xe=xk*(eh%*%matrix(1,1,k)) # x*e # omega=t(xe)%*%xe # Central part of Variance matrix # v=xxi%*%omega%*%xxi # variance matrix estimator # se=sqrt(diag(v)) # Standard errors # betas=cbind(beta,se) # estimates + standard errors # yf=xf%*%beta # Point forecast = Estimated change # print("Last 4 values, level and change") print(cbind(r[(n-2):(n+1)],dr[(n-3):n])) print("AR(12) Model") print("Coefficient Estimates and Standard Errors") print(betas) print("Point Forecast: Level and Change") print(cbind(r[n+1]+yf,yf)) # Model Selection # kn=kk+1 # Number of models = AR(0) through AR(kk) # yf=matrix(nrow=kn,ncol=1) # vector of forecasts (empty for now) # crit=matrix(nrow=kn,ncol=7) # vector of CV criterion (empty for now) # # Estimate Full Model for Mallows Criteria # xk=x[1:nn,] xxi=solve(t(xk)%*%xk) e=y-xk%*%xxi%*%(t(xk)%*%y) h=rowSums((xk%*%xxi)*xk) ek=e/(1-h) sigk=(t(e)%*%e)/nn 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) sig=(t(e)%*%e)/nn yf[k]=xf%*%beta aic=nn*log(sig)+k*2 # AIC # xe=xk*(eh%*%matrix(1,1,k)) # x*e # omega=t(xe)%*%xe # Central part of Variance matrix # b=matrix(1,1,k)%*%diag(xxi%*%omega) # correction for Robust Mallows # crit[k,1]=nn*log(sig)+k*log(nn) # BIC # crit[k,2]=aic # AIC # crit[k,3]=aic+2*(k+1)*(k+2)/(nn-k-2) # AICc # crit[k,4]=sig+2*k+sigk/nn # Mallows # crit[k,5]=sig+2*b/nn # Robust Mallows # crit[k,6]=sig+2*k+sig/nn # FPE # crit[k,7]=(t(eh)%*%eh)/nn # Cross-Validation # } k=which.min(crit[,7]) s=cbind(seq(0,kk),yf,crit[,1],crit[,3],crit[,5],crit[,7]) print("AR(0) through AR(24) Models, Point Forecasts, BIC, AICc, Robust Mallow, CV") print(s) print("Selected AR Order") print(k-1) print("Point Forecast: Level and Change") print(cbind(r[n+1]+yf[k],yf[k]))