########################################################################## # TSA2 ########################################################################## library(xts) library(fGarch) install.packages("devtools") install.packages("astsa") library(astsa) #################################################################################################### # Project with M1 Money Stock data, Weekly, Seasonally Adjusted from https://research.stlouisfed.org #################################################################################################### fed <- read.csv("/Users/wintenberger/DropSU/Lecture/LPSM/TSA/Test1/Data.csv", sep=";", dec=",") dim(fed)[1]->T Data<-ts(fed[,2],frequency=52,start = c(1975,2), end = c(2015,18)) ts.plot(Data) ########################################################################## # Stationarization: Log-ratios ########################################################################## data<-log(Data) ts.plot(data) data<-100*diff(log(Data)) ts.plot(data) anomalies<-which(abs(data)>6) data[anomalies]<-mean(data[anomalies+c(-1,1)]) ts.plot(data) test <- data[(length(data)-9):length(data)] data <- data[1:(length(data)-10)] ########################################################################## # Model Selection ########################################################################## p.opt = 10 q.opt = 13 fitarma<-arima(data,order=c(p.opt,0,q.opt),include.mean = FALSE) tsdiag(fitarma) ########################################################################### # Residuals as white noise but potential dependence in the squares ########################################################################### eps<-ts(fitarma$residuals,frequency=52,start = c(1975,2)) par(mfrow=c(1,1)) plot.ts(eps) acf(eps,lag.max=100) qqnorm(eps) qqline(eps,col=2) par(mfrow=c(1,2)) acf(eps) acf(eps^2) ########################################################################## # GARCH(1,1) estimation ########################################################################## ################## Qlik computation ##################################### objf.garch <- function(vartheta,eps,n,sig2init,r0=10){ omega <- vartheta[1] alpha <- vartheta[2] beta <- vartheta[3] sig2 <-rep(0,n) sig2[1] <- sig2init for(t in 2:n){ sig2[t] <- omega + alpha*eps[t-1]^2 + beta*sig2[t-1] } qml <- mean(eps[(r0+1):n]^2/sig2[(r0+1):n]+log(sig2[(r0+1):n])) qml } ################# Asymptotic variance approximation #################### VarAsymp<- function(omega,alpha,beta,eps,sig2init,petit,r0=10){ n <- length(eps) dersigma2<-matrix(0,nrow=3,ncol=n) sig2<-rep(0,n) sig2[1]<-sig2init for(t in 2:n){ vec<-c(1,eps[t-1]^2,sig2[t-1]) sig2[t]<-omega+beta*sig2[t-1]+alpha*eps[t-1]^2 dersigma2[1:3,t]<-vec/sig2[t]+beta*dersigma2[1:3,(t-1)] } eta <- eps[(r0+1):n]/sqrt(sig2)[(r0+1):n] eta <- eta/sd(eta) J<-dersigma2[1:3,(r0+1):n]%*%t(dersigma2[1:3,(r0+1):n])/(n-r0) kappa4<-mean(eta^4) {if(kappa(J)<1/petit) inv<-solve(J) else inv<-matrix(0,nrow=3,ncol=3)} var<-(kappa4-1)*inv list(var=var,residus=eta) } ######################### QMLE ################################################ estimGARCH<- function(omega,alpha,beta,eps,petit=sqrt(.Machine$double.eps),r0=10) { valinit<-c(omega,alpha,beta) n <- length(eps) sig2init<-var(eps[1:min(n,r0)]) res <- nlminb(valinit,objf.garch,lower=c(petit,0,0), upper=c(Inf,Inf,1), eps=eps,n=n,sig2init=sig2init) omega <- res$par[1] alpha<- res$par[2] beta <- res$par[3] var<-VarAsymp(omega,alpha,beta,eps,sig2init,petit=sqrt(.Machine$double.eps),r0=10) list(coef=c(omega,alpha,beta),residus=var$residus,var=var$var) } ################ Estimation ######################## par(mfrow=c(1,1)) omega.init<- 0.01 alpha.init<-0.01 beta.init<- 0.01 fitgarch<-estimGARCH(omega.init,alpha.init,beta.init,eps) par<-fitgarch$coef ################### Goodnees of fit diagnostic #################### res<-fitgarch$residus qqnorm(res) qqline(res,col=2) acf(res^2) ##################################################################### ################## Simple test procedure ########################### #################################################################### ##################### Test beta=0 ################################## beta=par[3] sigmabeta=fitgarch$var[3,3] se=sqrt(sigmabeta)/sqrt(T) t.value=beta/se p.value=pnorm(-t.value) # Symmetry print(p.value) ###################### One can reject beta=0 at a high level ####### ##################### For testing alpha = 0 one has to assume beta = 0 (ARCH) ######### ####################### Estimation ARCH(1) ############# estimARCH<- function(omega,alpha,beta=0,eps,petit=sqrt(.Machine$double.eps),r0=10) { valinit<-c(omega,alpha,beta) n <- length(eps) sig2init<-var(eps[1:min(n,5)]) res <- nlminb(valinit,objf.garch,lower=c(petit,0,0), upper=c(Inf,Inf,0), eps=eps,n=n,sig2init=sig2init) omega <- res$par[1] alpha<- res$par[2] var<-VarAsymp(omega,alpha,beta,eps,sig2init,petit=sqrt(.Machine$double.eps),r0=10) list(coef=c(omega,alpha),residus=var$residus,var=var$var[1:2,1:2]) } omega.init<- 0.01 alpha.init<-0.01 fitarch<-estimARCH(omega.init,alpha.init,eps=eps) par<-fitarch$coef res<-fitarch$residus qqnorm(res) qqline(res,col=2) acf(res^2) ##################### Test alpha=0 ###### alpha=par[2] sigmaalpha=fitarch$var[2,2] se=sqrt(sigmaalpha)/sqrt(T) t.value=alpha/se p.value=pnorm(-t.value) print(p.value) ##################### One can reject alpha = 0 at a high level ######## ######################################################################## #################### Prediction intervals ########################## ######################################################################## ##################### Two-step procedure ################## fitgarch<-estimGARCH(omega.init,alpha.init,beta.init,eps) sigmaGARCH <- fitgarch$coef[1]+fitgarch$coef[2]*eps[length(eps)]^2 + fitgarch$coef[3]*(eps[length(eps)]/fitgarch$residus[length(fitgarch$residus)])^2 sigmaGARCH <- rep(sigmaGARCH,10) for (i in 1:9){ sigmaGARCH[i+1] <- fitgarch$coef[1]+(fitgarch$coef[2]+fitgarch$coef[3])*sigmaGARCH[i] } se <- sqrt(predict(fitarma,n.ahead=10)$se^2 -fitarma$sigma2+sigmaGARCH) matplot(cbind(test,c(predict(fitarma,n.ahead=10)$pred),c(predict(fitarma,n.ahead=10)$pred+2*se),c(predict(fitarma,n.ahead=10)$pred-2*se)),type="l",lty=c(1,2,2,2),lwd=3,xlab="Horizon",ylab="Interval of prediction",main="Two-step procedure") ##################### One-step procedure ################# library(rugarch) model<-ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)), mean.model = list(armaOrder = c(p.opt, q.opt), include.mean = FALSE), distribution.model = "norm") modelfit<-ugarchfit(spec=model,data=data) par(mfrow=c(1,1)) plot(modelfit) forecast=ugarchforecast(modelfit,n.ahead = 10)@forecast matplot(cbind(test,c(forecast$series),c(forecast$series+2*forecast$sigma),c(forecast$series-2*forecast$sigma)),type="l",lty=c(1,2,2,2),lwd=3,xlab="Horizon",ylab="Interval of prediction",main="One-step procedure")