Simulating Risk-Neutral GARCH Paths for Option Pricing

# GARCH process simulation in R


# 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.

## Answer by Bob Jansen (score 1)

https://quant.stackexchange.com/a/75190

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.