################################################### # These functions are for estimating a GARCH(p,q) time series model in R using ranks # Code written by Huanhuan Wang, Department of Statistics, Northwestern University (HuanhuanWang2008@u.northwestern.edu) # supervised by Beth Andrews, Department of Statistics, Northwestern University (bandrews@northwestern.edu) # Requirement: the package "TSA" needs to be installed first install.packages("TSA") library ("TSA") ################################################### # Let "X" represent the observed GARCH series of length "n" # Here we will use a simulated GARCH(1,1) series of length n=1000 with parameters alpha0=0.1, alpha1=0.5 and beta1=0.4 # and iid N(0,1) noise n<-1000 X<-garch.sim(alpha=c(.1,.5), beta=c(.4), n = n, rnd = rnorm, ntrans = 100) ts.plot(X) #------------------------------------------------- # Choose GARCH model orders #-------------------------------------------------- p<-1 q<-1 ################################################### #--------------------------------------------------------------------------------- #Auxiliary functions: #--------------------------------------------------------------------------------- ################################################### #-------------------------------------------------- # the weight function lambda. # Here we have a weight function which is optimal when the GARCH noise is rescaled student's t with df=7 #------------------------------------------------- lambda <- function(x){ if(x<1) {c<-(7*(sqrt(5/7)*qt((x+1)/2,7))^2-5)/((sqrt(5/7)*qt((x+1)/2,7))^2+5)} if(x==1) {c<-7} return(c) } ################################################### ################################################### #-------------------------------------------------- # the weight function lambda squared. #------------------------------------------------- lambda2 <- function(x){ c <- lambda(x)^2 return(c) } ################################################### ################################################### #-------------------------------------------------- # getD returns the value of D for a candidate parameter vector theta #------------------------------------------------- getD <- function(theta){ #-------------------------------------------------- # calculate sigma square sigmasq <- array(1, (n+max(0,q-p))) for (t in (p+1):n){ sigmasq[t+max(0,q-p)]<-1+sum(theta[1:p]*((X[(t-1):(t-p)])^2)) if(q>0){sigmasq[t+max(0,q-p)]<-sigmasq[t+max(0,q-p)]+sum(theta[(p+1):(p+q)]*sigmasq[(t+max(0,q-p)-1):(t+max(0,q-p)-q)])} } #-------------------------------------------------- # calculate epsilon epsilon <- log(X^2)-log(abs(sigmasq[(1+max(0,q-p)):(n+max(0,q-p))])) #-------------------------------------------------- # sort epsilon sortep <- sort(epsilon[p+1:n]) #-------------------------------------------------- # calculate D D<-0 for (t in 1:(n-p)){ D<-D+lambda(t/(n-p+1))*(sortep[t]-mean(sortep)) } if(min(theta)<0){D=D+10000000} if((q>0) & (sum(theta[(p+1):(p+q)])>=1)){D=D+10000000} return(D) } ################################################### #--------------------------------------------------------------------------------- #End of auxiliary functions #--------------------------------------------------------------------------------- ################################################### #-------------------------------------------------- # Main Body #-------------------------------------------------- # generate 100 initial values for theta=(alpha1/alpha0, ... , alphap/alpha0, beta1, ... ,betaq) # under the condition that the sum of alpha1, ... ,alphap, beta1, ... ,betaq is less than 1 #-------------------------------------------------- no_initial <- 100 temp_alphabeta <-matrix(0,nr=no_initial,nc=p+q+1) temp_theta <-matrix(0,nr=no_initial,nc=p+q) count<-0 repeat{ alphabeta <- runif(p+q) if (sum(alphabeta)<1){ count<-count + 1 temp_alphabeta[count,] <- c(var(X)*(1-sum(alphabeta)),alphabeta) } if(count == no_initial)break } for (i in 1:p){ temp_theta[,i]<-temp_alphabeta[,(i+1)]/temp_alphabeta[,1] } if (q>0){ for (i in (p+1):(p+q)){ temp_theta[,i]<-temp_alphabeta[,(i+1)] } } #-------------------------------------------------- # calculate D for the initial parameter vectors theta #-------------------------------------------------- temp_D<- array(1000,no_initial) for (i in 1:no_initial){ temp_D[i] <-getD(temp_theta[i,]) } #-------------------------------------------------- # find the 3 initial values with the smallest values for D #-------------------------------------------------- store_D <- array(1000,3) store_theta <-matrix(0,nr=3,nc=p+q) for (j in 1:3){ store_D[j] <- min(temp_D) store_theta[j,] <- temp_theta[which.min(temp_D),] temp_D[which.min(temp_D)] <- max(temp_D)+1 } #-------------------------------------------------- # Using the 3 initial values as starting points, find optimized values for theta # The Rank estimate has the smallest corresponding D value #-------------------------------------------------- min_theta <- store_theta[1,] min_D <- store_D[1] for (j in 1:3){ theta<-store_theta[j,] theta<-optim(par=theta,fn=getD) tempD<-theta\$value if (tempD0){sigma_sq[t+max(0,q-p)]<-sigma_sq[t+max(0,q-p)]+sum(min_theta[(p+1):(p+q)]*sigma_sq[(t+max(0,q-p)-1):(t+max(0,q-p)-q)])} } ep<-log((X[(p+1):n])^2)-log(sigma_sq[(p+1+max(0,q-p)):(n+max(0,q-p))]) sortep<-sort(ep) b<-.9*n^(-.2)*min(sd(ep),(quantile(ep,.75)-quantile(ep,.25))/1.34) K<-0 for(t in 1:(n-p)){ for(s in 1:(n-p)){ K<-K+dnorm((ep[s]-sortep[t])/b)*(lambda(t/(n-p))-lambda((t-1)/(n-p)))/(b*n) } } #-------------------------------------------------- # Estimate Gamma #-------------------------------------------------- sigma_sq_der<-matrix(0,nr=n+max(0,q-p),nc=p+q) for(t in (p+1):n) { if(q==0) { for(i in 1:p) { sigma_sq_der[t,i]=X[t-i]^2 } } if(q>0) { for(i in 1:p) { sigma_sq_der[t+max(0,q-p),i]=X[t-i]^2+sum(min_theta[(p+1):(p+q)]*sigma_sq_der[(t+max(0,q-p)-1):(t+max(0,q-p)-q),i]) } for(i in (p+1):(p+q)){ sigma_sq_der[t+max(0,q-p),i]=sigma_sq[t+max(0,q-p)-i+p]+sum(min_theta[(p+1):(p+q)]*sigma_sq_der[(t+max(0,q-p)-1):(t+max(0,q-p)-q),i]) } } } ep_der<-matrix(0,nr=n+max(0,q-p),nc=p+q) for(i in 1:(p+q)) { ep_der[,i]=-sigma_sq_der[,i]/sigma_sq } Gamma=matrix(0,nr=p+q,nc=p+q) for(i in 1:(p+q)) { for(j in 1:(p+q)) { m1=sum(ep_der[(p+1+max(0,q-p)):(n+max(0,q-p)),i])/n m2=sum(ep_der[(p+1+max(0,q-p)):(n+max(0,q-p)),j])/n Gamma[i,j]=sum((ep_der[(p+1+max(0,q-p)):(n+max(0,q-p)),i]-m1)*(ep_der[(p+1+max(0,q-p)):(n+max(0,q-p)),j]-m2))/n } } #-------------------------------------------------- # Estimate Sigma #-------------------------------------------------- Sigma<-J*K^(-2)*solve(Gamma) #-------------------------------------------------- # Compute Wald statistics for testing H_0:theta_i=0, assuming all other parameters are positive # and display results #-------------------------------------------------- Wald<-array(0,p+q) for(i in 1:(p+q)) { Wald[i]<-n*min_theta[i]^2/Sigma[i,i] } cat ("Wald statistics for testing H_0:theta_i=0, assuming all other parameters are positive:",Wald,"\n") cat ("Reject H_0:theta_i=0 at .05 level of significance if test stat > ",qchisq(1-2*.05,df=1),"\n") #-------------------------------------------------- #-------------------------------------------------- # Assuming all parameters are positive, display 95% confidence intervals for the elements of theta #-------------------------------------------------- cat("Assuming all parameters are positive, 95% confidence intervals for the elements of theta:","\n") for(i in 1:(p+q)) { left<-min_theta[i]-1.96*sqrt(Sigma[i,i]/n) right<-min_theta[i]+1.96*sqrt(Sigma[i,i]/n) cat ("point estimate:",min_theta[i],", 95% CI: (",max(0,left),",",right,")","\n") } #-------------------------------------------------- ################################################### # Residual Analysis #-------------------------------------------------- #-------------------------------------------------- # GARCH residuals ("Z hat") standardized to have variance one: #-------------------------------------------------- sigma_sq <- array(1, n+max(0,q-p)) for (t in (p+1):n){ sigma_sq[t+max(0,q-p)]<-1+sum(min_theta[1:p]*((X[(t-1):(t-p)])^2)) if(q>0){sigma_sq[t+max(0,q-p)]<-sigma_sq[t+max(0,q-p)]+sum(min_theta[(p+1):(p+q)]*sigma_sq[(t+max(0,q-p)-1):(t+max(0,q-p)-q)])} } Z<-X/sqrt(sigma_sq[(1+max(0,q-p)):(n+max(0,q-p))]) Z<-Z/sd(Z) ts.plot(Z) #-------------------------------------------------- #Sample ACFs of residuals, absolute values of residuals, and squared residuals #-------------------------------------------------- acf(Z) acf(abs(Z)) acf(Z^2) #-------------------------------------------------- #Consider shape of noise distribution #-------------------------------------------------- qqnorm(Z) #Normal qq-plot of residuals hist(Z) #histogram of residuals plot(density(Z)) #kernel density estimate #-------------------------------------------------- ###################################################