|
| 1 | +#get the initial value of beta |
| 2 | +datamatrix=cbind(DBframe$index,DBframe$PE,DBframe$turnover,DBframe$csv) |
| 3 | +colnames(datamatrix)=c("index","PE","turnover","csv") |
| 4 | + |
| 5 | + |
| 6 | +Y=DBframe$Y |
| 7 | + |
| 8 | +result=lm(formula = Y ~ datamatrix) |
| 9 | +summary(result) |
| 10 | + |
| 11 | +#test heteroscedasticity |
| 12 | +par(mfrow=c(2,2)) # init 4 charts in 1 panel |
| 13 | +plot(result) |
| 14 | +library("lmtest") |
| 15 | +restest=bptest(result) |
| 16 | + |
| 17 | +#get the initial beta |
| 18 | +beta0=result$coefficients |
| 19 | + |
| 20 | +#estimate alpha in ARCH |
| 21 | +res = result$resid |
| 22 | +res2=res^2 |
| 23 | +res=as.matrix(res2) |
| 24 | + |
| 25 | +m=length(res) |
| 26 | +U=res[4:m] |
| 27 | +v1=matrix(res[3:(m-1)],nrow=m-3) |
| 28 | +v2=matrix(res[2:(m-2)],nrow=m-3) |
| 29 | +v3=matrix(res[1:(m-3)],nrow=m-3) |
| 30 | + |
| 31 | +V=cbind(v1,v2,v3) |
| 32 | + |
| 33 | +lm2=lm(U~V) |
| 34 | +summary(lm2) |
| 35 | + |
| 36 | +alphat=lm2$coefficients |
| 37 | + |
| 38 | +#calculate sigma |
| 39 | +v4=matrix(rep(1,m-3),nrow=m-3) |
| 40 | + |
| 41 | +Vt=cbind(v4,v1,v2,v3) |
| 42 | + |
| 43 | +sigma=Vt%*%alphat |
| 44 | + |
| 45 | +#calculate omega |
| 46 | +et=result$resid |
| 47 | +et=as.vector(et) |
| 48 | + |
| 49 | + |
| 50 | +et2=matrix(0,nrow=3896) |
| 51 | +mode(et2) |
| 52 | +et2[which(et<0)]=1 |
| 53 | +et2[which(et==0)]=1 |
| 54 | +et2[which(et>0)]=0 |
| 55 | +et2=as.vector(et2) |
| 56 | + |
| 57 | + |
| 58 | +thetavec=rep(theta,3896) |
| 59 | +thetavec=as.vector(thetavec) |
| 60 | + |
| 61 | +omega=abs(thetavec-et2) |
| 62 | +is.vector(omega) |
| 63 | +omega=omega[4:3896] |
| 64 | + |
| 65 | + |
| 66 | +#calculate the new beta |
| 67 | +Xt=cbind(1,datamatrix[4:3896,]) |
| 68 | +j=1 |
| 69 | +xxsummatrix=matrix(0,nrow=5,ncol = 5) |
| 70 | +xysummatrix=matrix(0,nrow = 5,ncol = 1) |
| 71 | +ynew=Y[4:3896] |
| 72 | +for (j in 1:3893) { |
| 73 | + xxsummatrix=xxsummatrix+omega[j]/sigma[j]*crossprod(t(Xt[j,]),Xt[j,]) |
| 74 | + xysummatrix=xysummatrix+omega[j]/sigma[j]*crossprod(t(Xt[j,]),ynew[j]) |
| 75 | + } |
| 76 | + |
| 77 | +solve(xxsummatrix) |
| 78 | +beta1=crossprod(t(solve(xxsummatrix)),xysummatrix) |
| 79 | + |
| 80 | + |
| 81 | + |
0 commit comments