Volatility updating with Heston-Nandi model 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} &amp;=&amp; r + \gamma h_{t+1} + \sqrt{h_{t+1}} z_{t+1}\\
    h_{t+1} &amp;=&amp; \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 &gt; 0,\ \beta_1 \ge 0,\ \beta_2 \ge 0\$, and \$\gamma &gt; -\frac{1}{2}\$.


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


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

$$
\left\{
  \begin{array}{rcl}
    h_{t+1} &amp;=&amp; \beta_0 + \beta_1 h_t + \beta_2\left(z_t - \beta_3 \sqrt{h_t}\right)^2 \\
    z_t     &amp;=&amp; \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&lt;-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 &lt;- 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?
First, in R we use &lt;- 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 &lt;-  rep(P_0, T) you could have P &lt;- vector(mode = &quot;numeric&quot;, length = T).
This is just a suggestion: But there is no reason to not name your parameters:
para &lt;- c(0.1, 0.2, 0.3, 0.4, 0.5, 0.7) # (beta_0,beta_1, beta_2, beta_3, r, gamma)
para &lt;- 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[&quot;gamma&quot;], but not para$gamma.
Also, there is no need for para_vol&lt;-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 &lt;-  3000
# For the example I simulate a GARCH
#process parameters
eta &lt;-  0.2 #eta = 0 is equivalent to Geometric Brownian Motion
mu &lt;-  100 #the mean of the process

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

# library(ggfortify)
# autoplot(sigma)

P_0 &lt;- mu #starting price, known
# P &lt;-  rep(P_0, T)
P &lt;- vector(mode = &quot;numeric&quot;, length = T)
P[1] &lt;- P_0

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

# Set the parameters :
para &lt;-
  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 &lt;- vector(&quot;numeric&quot;, length = length(P))
vol[1] &lt;-  sd(P)
for (i in 2:length(P)) {
  # para_vol &lt;- para[1:6]
  vol[i] &lt;-
    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)