GARCH(1,1) with Student’s \(t\) errors

Hamiltonian MC via STAN

1 Downloading stock returns data

#install.packages("rstan")
#install.packages("quantmod")
options(mc.cores=4)
options(stringsAsFactors = FALSE)
library("rstan")

## Loading required package: StanHeaders

## 
## rstan version 2.32.7 (Stan version 2.32.2)

## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)

library("quantmod")

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## Loading required package: TTR

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

getSymbols("AMZN", src = "yahoo", from = "2019-01-01", to = "2025-05-16")

## [1] "AMZN"

amazon = as.double(AMZN$AMZN.Adjusted)
getSymbols("PBR", src = "yahoo", from = "2019-01-01", to = "2025-05-16")

## [1] "PBR"

pbr = as.double(PBR$PBR.Adjusted)
getSymbols("^GSPC", src = "yahoo", from = "2019-01-01", to = "2025-05-16")

## [1] "GSPC"

sp500 = as.double(GSPC$GSPC.Adjusted)
getSymbols("^BVSP", src = "yahoo", from = "2019-01-01", to = "2025-05-16")

## [1] "BVSP"

ibovespa = as.double(BVSP$BVSP.Adjusted)

1.1 Closing prices

par(mfrow=c(1,1))
ts.plot(amazon/amazon[1],ylim=c(0.3,3.5))
lines(pbr/pbr[1],col=2)
lines(sp500/sp500[1],col=3)
lines(ibovespa/ibovespa[1],col=4)
abline(h=1,lty=2)
legend("topleft",legend=c("Amazon","Petrobrás","SP500","Ibovespa"),col=1:4,lwd=2,lty=1,bty="n")

1.2 Returns

par(mfrow=c(2,2))
ts.plot(diff(log(amazon)),ylim=c(-0.16,0.16),main="Amazon",ylab="Returns")
abline(h=0,col=2)
ts.plot(diff(log(pbr)),main="Petrobrás",ylab="Returns")
abline(h=-0.15,lty=2)
abline(h=0.15,lty=2)
abline(h=0,col=2)
ts.plot(diff(log(sp500)),ylim=c(-0.16,0.16),main="SP500",ylab="Returns")
abline(h=0,col=2)
ts.plot(diff(log(ibovespa)),ylim=c(-0.16,0.16),main="Ibovespa",ylab="Returns")
abline(h=0,col=2)

2 Petrobrás

n        = length(pbr)
ret      = pbr[2:n]/pbr[1:(n-1)]-1

par(mfrow=c(1,1))
acf(ret,lag=100,main="Log returns",ylim=c(-0.2,0.2))

acf(ret^2,lag=100,main="Squared log returns",ylim=c(-0.1,0.4))

2.1 Stan code

# Filename: garch.stan
data {
  int<lower=0> n;
  real y[n];
}

parameters {
  real<lower=0> alpha0;
  real<lower=0,upper=1> alpha1;
  real<lower=0,upper=(1-alpha1)> beta1;
  real<lower=0> sigma1;
  real<lower=2> nu;
}

transformed parameters {
  real<lower=0> sigma[n];
  sigma[1] = sigma1;
  for (t in 2:n)
    sigma[t]=sqrt(alpha0+alpha1*pow(y[t-1],2)+beta1*pow(sigma[t-1],2));
}

model {
  y ~ student_t(nu,0.0,sigma);
}

2.2 Running stan code for GARCH(1,1)

y      = ret
y      = (y-mean(y))/sqrt(var(y))
n      = length(y)
model  = stan_model("garch.stan")
fit    = sampling(model,data=list(n=n,y=y),iter=500,chains=4)
params = extract(fit)
qsig   = t(apply(params$sigma,2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(2,4))
ts.plot(params$alpha0,ylab="",xlab="Iteration",main=expression(alpha[0]))
ts.plot(params$alpha1,ylab="",xlab="Iteration",main=expression(alpha[1]))
ts.plot(params$beta1,ylab="",xlab="Iteration",main=expression(beta[1]))
ts.plot(params$nu,ylab="",xlab="Iteration",main=expression(nu))
hist(params$alpha0,prob=TRUE,main="",xlab="")
hist(params$alpha1,prob=TRUE,main="",xlab="")
hist(params$beta1,prob=TRUE,main="",xlab="")
hist(params$nu,prob=TRUE,main="",xlab="")

2.3 Posterior of \(\alpha_1+\beta_1\)

par(mfrow=c(1,1))
hist(params$alpha1+params$beta1,prob=TRUE,xlab="",main=expression(alpha[1]+beta[1]))

2.4 Time-varying standard deviations

par(mfrow=c(1,1))
plot(abs(y),type="h",col=grey(0.8),xlab="Days",ylab="Log returns (absolute value)")
lines(qsig[,1],col=3,lwd=3)
lines(qsig[,3],col=3,lwd=3)
lines(qsig[,2],col=2)
abline(h=1,lty=1)

3 Comparison

price   = amazon
n       = length(price)
y       = price[2:n]/price[1:(n-1)]-1
y       = (y-mean(y))/sqrt(var(y))
n       = length(y)
model   = stan_model("GARCH.stan")
fit     = sampling(model,data=list(n=n,y=y),iter=500,chains=4)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

par1    = extract(fit)

price   = pbr
n       = length(price)
y       = price[2:n]/price[1:(n-1)]-1
y       = (y-mean(y))/sqrt(var(y))
n       = length(y)
model   = stan_model("GARCH.stan")
fit     = sampling(model,data=list(n=n,y=y),iter=500,chains=4)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

par2    = extract(fit)

price   = sp500
n       = length(price)
y       = price[2:n]/price[1:(n-1)]-1
y       = (y-mean(y))/sqrt(var(y))
n       = length(y)
model   = stan_model("GARCH.stan")
fit     = sampling(model,data=list(n=n,y=y),iter=500,chains=4)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

par3    = extract(fit)

price   = ibovespa
n       = length(price)
y       = price[2:n]/price[1:(n-1)]-1
y       = (y-mean(y))/sqrt(var(y))
n       = length(y)
model   = stan_model("GARCH.stan")
fit     = sampling(model,data=list(n=n,y=y),iter=500,chains=4)
par4    = extract(fit)

tab1 = apply(cbind(par1$alpha0,par1$alpha1,par1$beta1,par1$nu),2,quantile,c(0.05,0.5,0.95))
tab2 = apply(cbind(par2$alpha0,par2$alpha1,par2$beta1,par2$nu),2,quantile,c(0.05,0.5,0.95))
tab3 = apply(cbind(par3$alpha0,par3$alpha1,par3$beta1,par3$nu),2,quantile,c(0.05,0.5,0.95))
tab4 = apply(cbind(par4$alpha0,par4$alpha1,par4$beta1,par4$nu),2,quantile,c(0.05,0.5,0.95))

tab.alpha0 = cbind(tab1[,1],tab2[,1],tab3[,1],tab4[,1])
tab.alpha1 = cbind(tab1[,2],tab2[,2],tab3[,2],tab4[,2])
tab.beta1  = cbind(tab1[,3],tab2[,3],tab3[,3],tab4[,3])
tab.nu     = cbind(tab1[,4],tab2[,4],tab3[,4],tab4[,4])
colnames(tab.alpha0) = c("apple","petrobras","sp500","ibovespa")
colnames(tab.alpha1) = c("apple","petrobras","sp500","ibovespa")
colnames(tab.beta1)  = c("apple","petrobras","sp500","ibovespa")
colnames(tab.nu)     = c("apple","petrobras","sp500","ibovespa")

3.1 \(\alpha_0\)

tab.alpha0

##           apple   petrobras      sp500    ibovespa
## 5%  0.004352389 0.004676527 0.01095361 0.008952082
## 50% 0.013944802 0.011288371 0.01729598 0.017190574
## 95% 0.032134317 0.024786102 0.02632215 0.028877962

3.2 \(\alpha_1\)

tab.alpha1

##          apple  petrobras      sp500   ibovespa
## 5%  0.02880990 0.01985034 0.09586442 0.04744191
## 50% 0.05193612 0.03294889 0.12379141 0.06571933
## 95% 0.09113641 0.05393075 0.15666511 0.09124835

3.3 \(\beta_1\)

tab.beta1

##         apple petrobras     sp500  ibovespa
## 5%  0.8237144 0.8615268 0.7610315 0.8476439
## 50% 0.8994968 0.9192599 0.8060386 0.8902761
## 95% 0.9439138 0.9520891 0.8428494 0.9194046

3.4 \(\nu\)

tab.nu

##        apple petrobras    sp500  ibovespa
## 5%  4.398443  3.901891 5.460582  7.732055
## 50% 5.428589  4.693051 7.010412 10.697054
## 95% 6.890187  5.645143 9.634685 16.126130

3.5 \(\alpha_1+\beta_1\)

tab.alpha1+tab.beta1

##         apple petrobras     sp500  ibovespa
## 5%  0.8525243 0.8813771 0.8568959 0.8950858
## 50% 0.9514329 0.9522088 0.9298300 0.9559954
## 95% 1.0350502 1.0060198 0.9995145 1.0106530