forked from burakbayramli/books
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rates1.r
98 lines (85 loc) · 3.08 KB
/
rates1.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
# 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)]
X11()
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]))