# rates3.r # library(quadprog) # You may need to install the package # library(tseries) # You may need to install the package # library(quantreg) # 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=r[n+1]+t(yf)%*%w x.arch <- garch(e,order=c(1,1)) archc=coef(x.arch) sd=predict(x.arch) like=logLik(x.arch) sd <- as.matrix(sd[,1]) var <- as.matrix(sd^2) sdf=sqrt(archc[1]+archc[2]*(e[nn]^2)+archc[3]*var[nn,1]) # Quantile Error Calculation # ep=e/sd print(summary(ep)) q1=coef(rq(ep~1,.1)) q2=coef(rq(ep~1,.25)) q3=coef(rq(ep~1,.75)) q4=coef(rq(ep~1,.9)) i1=yff+sdf*q1 i2=yff+sdf*q2 i3=yff+sdf*q3 i4=yff+sdf*q4 print("10%, 25%, 75%, 90% quantiles of normalized residuals") print(cbind(q1,q2,q3,q4)) print("Point Forecast") print(yff) print("50% Forecast Interval") print(cbind(i2,i3)) print("80% Forecast Interval") print(cbind(i1,i4)) # Quantile Regression approach to forecast intervals # k=3 xk=x[1:nn,1:k] xf=x[nn+1,1:k] x2=xk[,2:k] beta1=coef(rq(y~x2,.1)) beta2=coef(rq(y~x2,.25)) beta3=coef(rq(y~x2,.75)) beta4=coef(rq(y~x2,.9)) yf1=r[n+1]+xf%*%beta1 yf2=r[n+1]+xf%*%beta2 yf3=r[n+1]+xf%*%beta3 yf4=r[n+1]+xf%*%beta4 print("Quantile Regression Intervals") print("50% Forecast Interval") print(cbind(yf2,yf3)) print("80% Forecast Interval") print(cbind(yf1,yf4)) print("Estimated Coefficients") print(beta1) print(beta2) print(beta3) print(beta4) # Forecast Distribution # e1=as.matrix(sort(ep)) # Sort the normalized residuals # n1=nrow(e1) yp=matrix(1,n1,1)%*%yff+e1%*%sdf # Forecast Distribution # edf=seq(1,n1)/n1 plot(yp,edf,main="Interest Rate Forecast Distribution",type="l",xlab="",ylab="")