# GARCH process simulation in R

I'm trying to learn how to simulate the GARCH(1,1) for option pricing using Monte Carlo. I need to learn how to code the equations for the stock log returns and the variance process. I'm trying to reproduce the simple example given in Duan (2000), where all the information is given even the randomness for reproduction purposes. However, I can't really get the exact answer as in the given example. I think that I'm not coding the equations correctly. Below is my code


S0 <- 51    #current stock price
r <- 0.05   # interest rate
sd1 <- 0.2  # initial standard deviation
K <- 50     # strike price
T2M <- 2    # 2 days to maturity
sim <- 10   # number of simulations

beta0 <- 1e-05
beta1 <- 0.8
beta2 <- 0.1
theta <- 0.5
lambda <- 0.3

# given randomness term for reproducible simulation
z1 <- c(-0.8131,-0.547,0.4109,0.437,0.5413,-1.0472,0.3697,-2.0435,-0.2428,0.3091)
z2 <- c(0.7647, 0.5537, 0.0835, -0.6313, -0.1772, 2.4048, 0.0706, -1.4961,-1.376,0.3845)

# error term
Z <- unname(cbind(z1,z2))

# error under risk neutral probability Q
ZQ <- Z + lambda

# one-step variance need previous step as input (GARCH(1,1)) under Q
variance <- function(beta0, beta1, beta2, theta, lambda, z, Sigma){
return(beta0 + beta1 * Sigma^2 + beta2*Sigma^2*(z -theta -lambda)^2)
}

# one-step stock process
S_T <- function(S0, r, Sigma, T2M, z){
# It returns an array of terminal stock prices.
return(S0*exp((r + lambda*Sigma - Sigma^2/2) + Sigma*sqrt(T2M)*z))
}

# payoff function for European call
payoff_call <- function(S_T, K, r, T2M){
# It returns an array which has the value of the call for each terminal stock price.
return(exp(-r*T2M)*pmax(S_T-K, 0) )
}

# calculate one-step stock price
S1 <- S_T(S0, r, sd1, T2M=T2M/2, z=ZQ[,1])
S1

# calculate standard deviation at step 2
sd2 <- sqrt(variance(beta0, beta1, beta2, theta, lambda, ZQ[,1], sd1))
sd2


My output for the standard deviation is

0.1972484 0.1907743 0.1790021 0.1789577 0.1789325 0.2039248 0.1791031 0.2405984 0.1849784 0.1793203


and stock at time 1 is

50.36042 53.11321 64.32868 64.66536 66.02844 48.05690 63.80079 39.37478 56.44494 63.03220


whereas in the example the author have different output

I would appreciate the help.

• It would be better if you have summarized in sufficient detail the analysis you are trying reproduce so that one can still understand your question should the link to the paper you have referred to is broken or changed. In addition, most people may not have time to go and read that paper to understand what you are trying to fully accomplish. Apr 11, 2023 at 1:13
• @Richard thank you for only making the image appear and not changing anything else from my original post. Just to clarify to others that Richard didn't add anything scientific to the original post. Apr 13, 2023 at 10:47

To be honest, I don't quite understand the steps in the linked PDF. The std. normals don't like standard normals at all, the absolute values are too high.

I did manage to match your stock prices with only two minor changes:

# error under risk neutral probability Q
ZQ <- Z - lambda # was plus


The interest rate and volatility are annualized but the simulations are based on days. Often, annualization is done based on 250 or 252 days, in this case 365 days are used. The code becomes

days_per_year <- 365
S1 <- S_T(S0, r / days_per_year, sd1 / sqrt(days_per_year), T2M=T2M/2, z=ZQ[,1])


This last issue is without a doubt a bug in the original code. 5% daily interest doesn't make sense.

• Thank you very much @Bob. lambda is the risk premium. It is the $(\mu-r)/\sigma$ in Example 2 Risk-Neutral Measure. According to this link, $Z^Q = Z + \lambda$ is the correct form. Therefore, it's not the right answer. Apr 12, 2023 at 20:39
• I think there is a confusion of terms somewhere. Your S_T function seems to correspond with the first equation in the linked PDF. This appears to be for simulation under the real world measure but you pass in ZQ and not Z. Under what measure are z1 and z2 simulated? Apr 13, 2023 at 5:54
• S_T should be under Q measure for simulation. However, in the paper I'm not sure about the assumption regarding z1 and z2. I assumed they're under P. Apr 13, 2023 at 10:42
• I think the simulation is under the real world measure as $dW_t$ is used and not $dW^*_t$. Apr 20, 2023 at 19:43