4
\$\begingroup\$

I'm trying to program a volatility updating rule using iteration. I start with the well-known Heston-Nandi model where the returns dynamics are:

$$ \left\{ \begin{array}{rcl} R_{t+1} &=& r + \gamma h_{t+1} + \sqrt{h_{t+1}} z_{t+1}\\ h_{t+1} &=& \beta_0 + \beta_1 h_t + \beta_2\left(z_t - \beta_3 \sqrt{h_t}\right)^2 \end{array} \right. $$

where

  • \$z_t\$ is an IID standard normal random variable
  • \$h_t\$ is time-varying squared volatility
  • \$\beta_0 > 0,\ \beta_1 \ge 0,\ \beta_2 \ge 0\$, and \$\gamma > -\frac{1}{2}\$.

I want to do is to write the code associate to the volatility updating rule, explained in this algorithm:

  1. Define \$h_0 = c\$ equals the given unconditional variance which is constant,
  2. Iteration for \$t \in \{1, 2, \cdots, n\}\$:

    $$ \left\{ \begin{array}{rcl} h_{t+1} &=& \beta_0 + \beta_1 h_t + \beta_2\left(z_t - \beta_3 \sqrt{h_t}\right)^2 \\ z_t &=& \frac{[R_t - r - \gamma h_t]}{\sqrt{h_t}} \end{array} \right. $$

    to obtain the returns based proxy for spot variances \$(h_t^R)_t\$. Which yields an updating function that exclusively involves observation:

    $$ h_{t+1} = \beta_0 + \beta_1 h_t + \frac{\beta_2}{h_t} \left( R_t - r - (\beta_3 + \gamma)h_t \right)^2$$

My program is the following:

library(fGarch)

T=3000
    # For the example I simulate a GARCH
    #process parameters
    eta = 0.2 #eta = 0 is equivalent to Geometric Brownian Motion
    mu = 100 #the mean of the process

    #GARCH volatility model
    specs = garchSpec(model = list(omega = 0.000001, alpha = 0.5, beta = 0.4)) 
    sigma = garchSim(spec = specs, n = T)

    P_0 = mu #starting price, known
    P = rep(P_0,T)

    for(i in 2:T){
      P[i] = P[i-1] + eta * (mu - P[i-1]) + sigma[i] * P[i-1]
    }

    # Set the parameters :
para<-c(0.1,0.2,0.3,0.4,0.5,0.7) # (beta_0,beta_1, beta_2, beta_3, r, gamma) 
    # Iteration to obtain the volatility associate to the model :

vol = c()
vol[1]=sd(P)
for (i in 2:length(P)){
      para_vol <- para[1:6]
      vol[i]=para_vol[1]+ (para_vol[2]*vol[i-1])+ (para_vol[3]/vol[i-1])*(P[i-1]-para_vol[5]-(para_vol[4]+para_vol[6])*vol[i-1])
      }
vol

This is an example where I simulate a GARCH (as data set). I am trying to extract the volatility associated with the Heston-Nandi model.

I know I'm using a lot of bad things for R, but I could not figure out a better solution. Are there any corrections or suggestions needed to improve this process?

\$\endgroup\$
1
1
\$\begingroup\$

First, in R we use <- for assignment. It has its advantages.

It is prefered that you don't use T and F as variable names as these are aliases to TRUE and FALSE.

Instead of this P <- rep(P_0, T) you could have P <- vector(mode = "numeric", length = T).

This is just a suggestion: But there is no reason to not name your parameters:

para <- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.7) # (beta_0,beta_1, beta_2, beta_3, r, gamma)
para <- c(beta_0 = 0.1,  beta_1 = 0.2, beta_2 = 0.3, beta_4 = 0.4, r = 0.5, gamma = 0.7) # (beta_0,beta_1, beta_2, beta_3, r, gamma)

Then, you may access them like para["gamma"], but not para$gamma. Also, there is no need for para_vol<-para[1:6] in the loop, but I wager this is remnant from other attempts.

In terms of performance, the pre-allocation I used did a 70 ms to 50 ms. Although, I'm not confident in those numbers, as I did it fairly unscientifically.

Also, the referenced post shows a more direct way of calculating this.

Collapsed version of the code I ended up with:

library(fGarch)

T <-  3000
# For the example I simulate a GARCH
#process parameters
eta <-  0.2 #eta = 0 is equivalent to Geometric Brownian Motion
mu <-  100 #the mean of the process

#GARCH volatility model
specs <-  garchSpec(model = list(
  omega = 0.000001,
  alpha = 0.5,
  beta = 0.4
))
sigma <-  garchSim(spec = specs, n = T)

# library(ggfortify)
# autoplot(sigma)

P_0 <- mu #starting price, known
# P <-  rep(P_0, T)
P <- vector(mode = "numeric", length = T)
P[1] <- P_0

for (i in 2:T) {
  P[i] <-  P[i - 1] + eta * (mu - P[i - 1]) + sigma[i] * P[i - 1]
}
# DOESN'T WORK
# P[-T] <-  P[-T] + eta * (mu - P[-T]) + sigma[-1] * P[-T]

# Set the parameters :
para <-
  c(0.1, 0.2, 0.3, 0.4, 0.5, 0.7) # (beta_0,beta_1, beta_2, beta_3, r, gamma)
# Iteration to obtain the volatility associate to the model :

vol <- vector("numeric", length = length(P))
vol[1] <-  sd(P)
for (i in 2:length(P)) {
  # para_vol <- para[1:6]
  vol[i] <-
    para_vol[1] + (para_vol[2] * vol[i - 1]) + (para_vol[3] /
                                                  vol[i - 1]) * (P[i - 1] - para_vol[5] - (para_vol[4] + para_vol[6]) * vol[i - 1])
}
vol

# autoplot(vol)
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.