Time-varying variance modeling

ARCH, GARCH, SV models

1 Obtaining the data

#install.packages("quantmod")
#install.packages("stochvol")
#install.packages("fGarch")
#install.packages("rugarch")
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

library(stochvol)
library("fGarch")

## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer
## attached to the search() path when 'fGarch' is attached.
## 
## If needed attach them yourself in your R script by e.g.,
##         require("timeSeries")

## 
## Attaching package: 'fGarch'

## The following object is masked from 'package:TTR':
## 
##     volatility

library("rugarch")

## Loading required package: parallel

# Define the ticker symbol for Petrobras Common ADRs on NYSE
ticker <- "PBR"

# Define the start and end dates
start_date <- "2000-01-01"
end_date <- Sys.Date() # Uses the current date

# Download the historical data
# The data is automatically assigned to an object named "PBR" (the ticker symbol)
getSymbols(
  Symbols = ticker,
  src = "yahoo",         # Source: Yahoo Finance is a common and reliable source
  from = start_date,
  to = end_date
)

## [1] "PBR"

# The result is an **xts** (eXtensible Time Series) object named PBR
head(PBR)

##            PBR.Open PBR.High  PBR.Low PBR.Close PBR.Volume PBR.Adjusted
## 2000-08-10 6.312500 7.156250 6.125000  7.156250  143929200     1.254157
## 2000-08-11 6.843750 7.281250 6.750000  7.140625   25139200     1.251419
## 2000-08-14 7.031250 7.218750 7.000000  7.109375    6935200     1.245943
## 2000-08-15 7.109375 7.359375 7.015625  7.140625    9078800     1.251419
## 2000-08-16 7.203125 7.500000 7.125000  7.328125   11728000     1.284280
## 2000-08-17 7.390625 7.781250 7.359375  7.687500    7148000     1.347261

# Check the structure and class of the data
str(PBR)

## An xts object on 2000-08-10 / 2025-11-28 containing: 
##   Data:    double [6364, 6]
##   Columns: PBR.Open, PBR.High, PBR.Low, PBR.Close, PBR.Volume ... with 1 more column
##   Index:   Date [6364] (TZ: "UTC")
##   xts Attributes:
##     $ src    : chr "yahoo"
##     $ updated: POSIXct[1:1], format: "2025-12-01 19:15:08"

class(PBR)

## [1] "xts" "zoo"

# View a summary of the data
summary(PBR)

##      Index               PBR.Open         PBR.High         PBR.Low      
##  Min.   :2000-08-10   Min.   : 2.388   Min.   : 2.575   Min.   : 2.335  
##  1st Qu.:2006-12-07   1st Qu.: 8.500   1st Qu.: 8.637   1st Qu.: 8.370  
##  Median :2013-04-08   Median :13.440   Median :13.650   Median :13.240  
##  Mean   :2013-04-05   Mean   :17.029   Mean   :17.293   Mean   :16.730  
##  3rd Qu.:2019-08-01   3rd Qu.:20.990   3rd Qu.:21.372   3rd Qu.:20.661  
##  Max.   :2025-11-28   Max.   :76.650   Max.   :77.610   Max.   :74.540  
##    PBR.Close        PBR.Volume         PBR.Adjusted    
##  Min.   : 2.388   Min.   :   296800   Min.   : 0.4556  
##  1st Qu.: 8.489   1st Qu.:  9435400   1st Qu.: 2.3789  
##  Median :13.455   Median : 15346000   Median : 4.4426  
##  Mean   :17.015   Mean   : 18115945   Mean   : 5.4132  
##  3rd Qu.:21.029   3rd Qu.: 23519825   3rd Qu.: 7.7709  
##  Max.   :75.190   Max.   :207564100   Max.   :18.2651

p = as.numeric(PBR$PBR.Adjusted)

# You can also generate a quick candlestick chart
# This demonstrates the power of the quantmod package
chartSeries(PBR$PBR.Adjusted,
            name = "Petrobras Prices",
            theme = "white",
            TA = NULL) # Remove default technical indicators for a cleaner view

chartSeries(diff(log(PBR$PBR.Adjusted)),
            name = "Petrobras Returns",
            theme = "white",
            TA = NULL) # Remove default technical indicators for a cleaner view

2 Statistical analysis

price = as.numeric(PBR$PBR.Adjusted)
n = length(price)
return = (price[2:n]/price[1:(n-1)]-1)
n = length(return)

par(mfrow=c(2,2))
ts.plot(price)
ts.plot(return)
acf(price)
acf(return^2)

3 ARCH model

fit.arch  = garchFit(~garch(1,0),data=return,trace=F,include.mean=FALSE)

fit.arch

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 0), data = return, include.mean = FALSE, 
##     trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 0)
## <environment: 0x11cbb8a38>
##  [data = return]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1  
## 0.00066851  0.33114049  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  6.685e-04   1.597e-05    41.87   <2e-16 ***
## alpha1 3.311e-01   2.345e-02    14.12   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  13382.38    normalized:  2.103155 
## 
## Description:
##  Mon Dec  1 19:15:09 2025 by user:

par(mfrow=c(1,1))
ts.plot(fit.arch@sigma.t,ylab="Standard deviation")
title("Petrobras - ARCH(1) modeling")

4 GARCH model

fit.garch = garchFit(~garch(1,1),data=return,trace=F,include.mean=F)    

fit.garch

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = return, include.mean = F, 
##     trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x11a5c9df8>
##  [data = return]
## 
## Conditional Distribution:
##  norm 
## 
## Coefficient(s):
##      omega      alpha1       beta1  
## 1.8817e-05  8.1406e-02  8.9841e-01  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  1.882e-05   2.963e-06    6.351 2.14e-10 ***
## alpha1 8.141e-02   7.065e-03   11.522  < 2e-16 ***
## beta1  8.984e-01   8.801e-03  102.080  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  13918.82    normalized:  2.187462 
## 
## Description:
##  Mon Dec  1 19:15:09 2025 by user:

par(mfrow=c(1,1))
yrange = range(fit.arch@sigma.t,fit.garch@sigma.t)
ts.plot(fit.arch@sigma.t,ylab="Standard deviation",ylim=yrange)
lines(fit.garch@sigma.t,col=2)
title("ARCH(1) vs GARCH(11)")

# GARCH with Student’s t errors

fit.garch.std = garchFit(~garch(1,1),cond.dist="std",data=return,trace=F,include.mean=F)

fit.garch.std

## 
## Title:
##  GARCH Modelling 
## 
## Call:
##  garchFit(formula = ~garch(1, 1), data = return, cond.dist = "std", 
##     include.mean = F, trace = F) 
## 
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## <environment: 0x11ab40858>
##  [data = return]
## 
## Conditional Distribution:
##  std 
## 
## Coefficient(s):
##      omega      alpha1       beta1       shape  
## 1.1404e-05  5.9208e-02  9.2773e-01  6.2342e+00  
## 
## Std. Errors:
##  based on Hessian 
## 
## Error Analysis:
##         Estimate  Std. Error  t value Pr(>|t|)    
## omega  1.140e-05   2.596e-06    4.392 1.12e-05 ***
## alpha1 5.921e-02   7.460e-03    7.937 2.00e-15 ***
## beta1  9.277e-01   9.052e-03  102.493  < 2e-16 ***
## shape  6.234e+00   4.568e-01   13.648  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log Likelihood:
##  14116.72    normalized:  2.218564 
## 
## Description:
##  Mon Dec  1 19:15:09 2025 by user:

par(mfrow=c(1,1))
yrange = range(fit.garch@sigma.t,fit.garch.std@sigma.t)
ts.plot(fit.garch@sigma.t,ylab="Standard deviation",ylim=yrange)
lines(fit.garch.std@sigma.t,col=2)
title("GARCH(1,1) vs GARCH(1,1)-Student t")

n0 = 1500

par(mfrow=c(1,1))
yrange = range(fit.garch@sigma.t[n0:(n-n0)],fit.garch.std@sigma.t[n0:(n-n0)])
plot(n0:(n-n0),fit.garch@sigma.t[n0:(n-n0)],ylab="Standard deviation",ylim=yrange,type="l",xlab="Days")
lines(n0:(n-n0),fit.garch.std@sigma.t[n0:(n-n0)],col=2)
title("GARCH(1,1) vs GARCH(1,1)-Student t")

5 SV model

ndraw = 10000
svfit = svsample(return,draws=ndraw)

## Argument 'y' (data vector) contains values very close to zero. I am applying an offset constant of size 3.12430360911304e-06 to do the auxiliary mixture sampling. If you want to avoid this, you might consider de-meaning the returns before calling this function.

## Done!

## Summarizing posterior draws...

param = svfit$para[[1]][,1:3]
stdev = exp(svfit$latent[[1]]/2)
qmcmc = t(apply(stdev,2,quantile,c(0.025,0.5,0.975)))

par(mfrow=c(2,2))
hist(param[,1],main="mu",prob=TRUE,xlab="")
hist(param[,2],main="phi",prob=TRUE,xlab="")
hist(param[,3],main="sigma",prob=TRUE,xlab="")
ts.plot(qmcmc[,2],ylab="Standard deviation")

par(mfrow=c(1,1))
ts.plot(qmcmc[,2],ylim=c(0,max(qmcmc[,2])),ylab="Standard deviation")
lines(0.1*abs(return),col=3,type="h")
legend("top",legend=c("Abs returns","tv-stdev"),col=c(3,1),lty=1,lwd=2,bty="n")
title("Petrobras - SV modeling")

par(mfrow=c(1,1))
plot((n-n0):n,qmcmc[(n-n0):n,2],ylim=c(-0.02,max(qmcmc[(n-n0):n,2])),
     ylab="Standard deviation",type="l",lwd=2,xlab="Day")
#lines((n-n0):n,0.2*abs(return[(n-n0):n]),col=2,type="h")
abline(h=sqrt(var(return[(n-n0):n])),lty=2,lwd=2,col=4)
abline(h=sqrt(var(return)),lty=4,lwd=2,col=4)
lines((n-n0):n,0.01+0.1*return[(n-n0):n],col=2)
legend("top",legend=c("returns","tv-stdev","stdev-1000","stdev-all"),
       col=c(2,1,4,4),lty=c(1,1,2,4),lwd=2,bty="n")
title("Petrobras - SV modeling")