rm(list=ls()) install.packages("fGarch") library(fGarch) library(stochvol) #====================================================== # Uploading SP500 data #====================================================== data = read.csv("sp500.csv",header=TRUE) n = nrow(data) price = data[,2] logprice = log(price) logreturn = diff(logprice) pdf(file="sp500-ts.pdf",width=6,height=9) par(mfrow=c(3,1)) ind = seq(1,n,length=5) plot(price,xlab="days",ylab="Price",main="",type="l",axes=FALSE) axis(2);axis(1,at=ind,lab=data[ind,1]);box() plot(logprice,xlab="days",ylab="Log-price",main="",type="l",axes=FALSE) axis(2);axis(1,at=ind,lab=data[ind,1]);box() ind = seq(2,n,length=5) plot(logreturn,xlab="days",ylab="Log-price",main="",type="l",axes=FALSE) axis(2);axis(1,at=ind,lab=data[ind,1]);box() dev.off() #====================================================== # ARCH(1) & GARCH(1,1) #====================================================== y = logreturn fit1 = garchFit(~garch(1,0),data=y,trace=F,include.mean=FALSE) fit2 = garchFit(~garch(1,1),data=y,trace=F,include.mean=FALSE) # Estimated coefficients fit1@fit$matcoef fit2@fit$matcoef # Residual analysis res1 = fit1@residuals/fit1@sigma.t res2 = fit2@residuals/fit2@sigma.t pdf(file="sp500-fit.pdf",width=9,height=6) par(mfrow=c(1,1)) range = range(fit1@sigma.t,fit2@sigma.t) plot(fit1@sigma.t,ylab="Standard deviation",xlab="Days",type="l",axes=FALSE,ylim=range) axis(2);axis(1,at=ind,lab=data[ind,1]);box() lines(fit2@sigma.t,col=2) legend("topleft",legend=c("ARCH(1)","GARCH(1,1)"),col=1:2,lty=1,lwd=2) dev.off() pdf(file="sp500-res.pdf",width=6,height=9) par(mfrow=c(3,2)) plot(res1,ylab="Standardized residuals",xlab="Days",type="l",axes=FALSE) axis(2);axis(1,at=ind,lab=data[ind,1]);box();title("ARCH(1)") plot(res2,ylab="Standardized residuals",xlab="Days",type="l",axes=FALSE) axis(2);axis(1,at=ind,lab=data[ind,1]);box();title("GARCH(1,1)") acf(res1^2,main="ACF squared residuals") acf(res2^2,main="ACF squared residuals") pacf(res1^2,main="PACF squared residuals") pacf(res2^2,main="PACF squared residuals") dev.off() Box.test(res1^2,lag=10,type='Ljung') Box.test(res2^2,lag=10,type='Ljung') #====================================================== # Stochastic Volatility AR(1) #====================================================== draws = svsample(y) apply(draws$para,2,summary) stdevs = t(apply(exp(draws$latent/2),2,quantile,c(0.025,0.5,0.975))) stdevm= apply(exp(draws$latent/2),2,mean) par(mfrow=c(1,1)) plot(fit2@sigma.t,ylab="Standard deviation",xlab="Days",type="l",axes=FALSE, ylim=c(0,max(fit2@sigma.t,stdevm))) axis(2);axis(1,at=ind,lab=data[ind,1]);box() lines(stdevm,col=2) lines(0.15*abs(y),type="h",col=grey(0.75)) legend("topleft",legend=c("SV-AR(1)","GARCH(1,1)","abs(y)"),col=c(2,1,grey(0.75)),bty="n",lty=1,lwd=2) #====================================================== # Unconditional variances and kurtosis #====================================================== alpha0 = fit1@fit$matcoef[1,1] alpha1 = fit1@fit$matcoef[2,1] uncond.var.arch = alpha0/(1-alpha1) uncond.kurt.arch= 3*((1-alpha1^2)/(1-3*alpha1^2)) alpha0 = fit2@fit$matcoef[1,1] alpha1 = fit2@fit$matcoef[2,1] beta1 = fit2@fit$matcoef[3,1] uncond.var.garch = alpha0/(1-alpha1-beta1) uncond.kurt.garch= 3*((1-(alpha1+beta1)^2)/(1-(alpha1+beta1)^2-2*alpha1^2)) uncond.var.sv = mean(exp(draws$para[,1])) uncond.kurt.sv = mean(3*exp(draws$para[,3]^2/(1-draws$para[,2]^2))) tab = data.frame( rbind(c(var(y),uncond.var.arch,uncond.var.garch,uncond.var.sv), sqrt(c(var(y),uncond.var.arch,uncond.var.garch,uncond.var.sv)), c(kurtosis(y)-3,uncond.kurt.arch,uncond.kurt.garch,uncond.kurt.sv))) tab = round(tab,7) colnames(tab) = c("Sample","ARCH(1)","GARCH(1,1)","SV-AR(1)") rownames(tab) = c("Variance","St.Dev.","Kurtosis") tab # GARCH(1,1) and SV-AR(1) for sub-groups # ================================ ini = c(1,2501,5001,7501,10001) fin = c(2500,5000,7500,10000,n-1) par.garch = matrix(0,5,3) par.sv = matrix(0,5,3) for (i in 1:5){ y1 = y[ini[i]:fin[i]] fit = garchFit(~garch(1,1),data=y1,trace=F,include.mean=FALSE) draws = svsample(y1) par.garch[i,] = fit@fit$par par.sv[i,] = apply(draws$para,2,mean) }