############################################################################### # # Bayesian GARCH(1,1) via Sampling Importance Resampling # # By Hedibert Freitas Lopes # # This version: May 9th 2017. # ############################################################################### # # Petroleo Brasileiro S.A. - Petrobras (PBR-A) # NYSE - NYSE Delayed Price. Currency in USD # https://finance.yahoo.com/quote/PBR-A # ############################################################################## rm(list=ls()) data = read.csv("PBR-A.csv",header=TRUE) n = nrow(data) ind = trunc(seq(1,n,length=6)) dates = data[ind,1] ind1 = trunc(seq(2,n,length=6)) dates1 = data[ind1,1] price = data[,6] return = 100*(data[2:n,6]/data[1:(n-1),6]-1) logreturn = log(return/100+1) pdf(file="petrobras-arch-garch.pdf",width=10,height=7) par(mfrow=c(1,1)) plot(price,xlab="Days",ylab="Price (US$)",main="",type="l",axes=FALSE) axis(2);box();axis(1,at=ind,lab=dates) par(mfrow=c(1,1)) acf(price,lag=740) par(mfrow=c(1,1)) plot(return,xlab="Days",ylab="Returns (%)",main="",type="l",axes=FALSE) axis(2);box();axis(1,at=ind1,lab=dates1) par(mfrow=c(1,1)) plot(logreturn,xlab="Days",ylab="Log returns",main="",type="l",axes=FALSE) axis(2);box();axis(1,at=ind1,lab=dates1) par(mfrow=c(1,2)) acf(logreturn,lag=100,main="Log returns") pacf(logreturn,lag=100,main="Log returns") # Are log returns Gaussian or heavier-tailed? ############################################# y = logreturn par(mfrow=c(1,1)) hist(y,xlab="Log returns",main="",prob=TRUE,breaks=seq(min(y),max(y),length=50)) box() den = density(y) x = density(y)$x yd = density(y)$y plot(x,yd,xlab="Losg returns",ylab="Density",type="l") lines(x,dnorm(x,mean(y),sqrt(var(y))),col=2) legend("topleft",legend=c("Log returns","Gaussian"),col=1:2,lty=2) plot(x,yd/dnorm(x,mean(y),sqrt(var(y))),xlim=c(-0.15,0.15),ylim=c(0,50),type="l", xlab="Log returns",ylab="Tail ratio") # Fitting ARCH(1) & GARCH(1,1) model to log-returns ############################################# par(mfrow=c(1,2)) acf(logreturn^2,lag=100,main="Squared log returns") pacf(logreturn^2,lag=100,main="Squared log returns") install.packages("fGarch") library(fGarch) fit.arch = garchFit(~garch(1,0),data=y,trace=F,include.mean=FALSE) fit.arch fit.arch@fit$matcoef par(mfrow=c(1,1)) plot(abs(y),xlab="Days",ylab="Absolute log returns",main="",type="h",axes=FALSE,col=gray(0.75)) axis(2);box();axis(1,at=ind1,lab=dates1) lines(fit.arch@sigma.t) fit.garch = garchFit(~garch(1,1),data=y,trace=F,include.mean=FALSE) fit.garch fit.garch@fit$matcoef par(mfrow=c(1,1)) plot(abs(y),xlab="Days",ylab="Absolute log returns",main="",type="h",axes=FALSE,col=gray(0.75)) axis(2);box();axis(1,at=ind1,lab=dates1) lines(fit.garch@sigma.t) loglike.garch11 = function(y,alpha0,alpha1,beta1){ n = length(y) h = alpha0 loglike = 0.0 for (t in 2:n){ h = alpha0+alpha1*y[t-1]^2+beta1*h loglike = loglike + dnorm(y[t],0,sqrt(h),log=TRUE) } return(loglike) } # Using the sampling importance resampling (SIR) Monte Carlo scheme # to sample from the joint posterior distribution of (alpha1,beta1) install.packages("mvtnorm") library("mvtnorm") alpha0 = 0.00001322076 alpha1 = 0.0786035 beta1 = 0.9086285 sd.alpha1 = 0.008726218 sd.beta1 = 0.009852661 rho = -0.75 mprior = c(alpha1,beta1) cov = rho*sd.alpha1*sd.beta1 vprior1 = matrix(c(sd.alpha1^2,cov,cov,sd.beta1^2),2,2) vprior = 10*vprior1 N = 30 alpha1s = seq(alpha1-4*sd.alpha1,alpha1+4*sd.alpha1,length=N) beta1s = seq(beta1-4*sd.beta1,beta1+4*sd.beta1,length=N) loglikes = matrix(0,N,N) logpriors = matrix(0,N,N) logprops = matrix(0,N,N) for (i in 1:N) for (j in 1:N){ loglikes[i,j] = loglike.garch11(y,alpha0,alpha1s[i],beta1s[j]) logpriors[i,j] = dmvnorm(c(alpha1s[i],beta1s[j]),mprior,vprior,log=TRUE) logprops[i,j] = dmvt(c(alpha1s[i],beta1s[j]),mprior,vprior1,df=5,log=TRUE) } logposts = loglikes+logpriors likes = exp(loglikes-max(loglikes)) priors = exp(logpriors-max(logpriors)) posts = exp(logposts-max(logposts)) props = exp(logprops-max(logprops)) par(mfrow=c(1,2)) contour(alpha1s,beta1s,likes,drawlabels=FALSE,xlab=expression(alpha[1]),ylab=expression(beta[1])) contour(alpha1s,beta1s,priors,col=2,add=TRUE,drawlabels=FALSE) title("Likelihood & prior contours") legend("topright",legend=c("Likelihood","Prior"),col=1:2,lty=1) contour(alpha1s,beta1s,posts,drawlabels=FALSE,xlab=expression(alpha[1]),ylab=expression(beta[1]),col=4) contour(alpha1s,beta1s,props,col=5,add=TRUE,drawlabels=FALSE) title("Posterior and proposal contours") legend("topright",legend=c("Posterior","Proposal"),col=c(4,5),lty=1) # SIR set.seed(2325) M = 1000 draws = matrix(mprior,M,2,byrow=TRUE)+rmvt(M,vprior1,df=5) w = rep(0,M) for (i in 1:M) w[i] = loglike.garch11(y,alpha0,draws[i,1],draws[i,2])+dmvnorm(c(draws[i,1],draws[i,2]),mprior,vprior,log=TRUE)-dmvt(c(draws[i,1],draws[i,2]),mprior,vprior1,df=5,log=TRUE) w = exp(w-max(w)) ii = sample(1:M,size=M,replace=TRUE,prob=w) draws1 = draws[ii,] par(mfrow=c(1,2)) contour(alpha1s,beta1s,props,drawlabels=FALSE,xlab=expression(alpha[1]),ylab=expression(beta[1])) points(draws,col=grey(0.75),pch=16) contour(alpha1s,beta1s,props,drawlabels=FALSE,add=TRUE) title("Proposal contours and draws") contour(alpha1s,beta1s,posts,drawlabels=FALSE,xlab=expression(alpha[1]),ylab=expression(beta[1])) points(draws1,col=grey(0.75),pch=16) contour(alpha1s,beta1s,posts,drawlabels=FALSE,add=TRUE) title("Posterior contours and draws") #dev.off() par(mfrow=c(1,2)) hist(draws1[,1],xlab="",prob=TRUE,main=expression(alpha[1])) hist(draws1[,2],xlab="",prob=TRUE,main=expression(beta[1])) par(mfrow=c(1,1)) hist(draws1[,1]+draws1[,2],xlab="",prob=TRUE,main=expression(alpha[1]+beta[1])) sig2 = matrix(0,M,n) for (t in 2:n) sig2[,t] = alpha0+draws1[,1]*y[t-1]^2+draws[,2]*sig2[,t-1] qsig2 = apply(sqrt(sig2),2,quantile,c(0.05,0.5,0.95)) par(mfrow=c(1,1)) plot(qsig2[2,2:n],xlab="Days",ylab="Log returns",main="",type="l",axes=FALSE,ylim=range(qsig2)) axis(2);box();axis(1,at=ind1,lab=dates1) lines(qsig2[1,2:n],col=2) lines(qsig2[3,2:n],col=2) par(mfrow=c(1,1)) plot(qsig2[2,(n-200):n],xlab="Last 200 days Days",ylab="Standard deviations",main="",ylim=range(qsig2[,(n-200):n]),pch=16) for (i in (n-200):n) segments(i-n+201,qsig2[1,i],i-n+201,qsig2[3,i],) title("July 20th 2016 to May 5th 2017" ) dev.off()