Simulación de procesos GARCH en R

Question

Estoy tratando de aprender a simular el GARCH(1,1) para la valoración de opciones utilizando Monte Carlo. Necesito aprender a codificar las ecuaciones para los retornos logarítmicos de las acciones y el proceso de varianza. Estoy tratando de reproducir el ejemplo simple dado en Duan (2000) donde se da toda la información, incluso la aleatoria, con fines de reproducción. Sin embargo, no puedo obtener la respuesta exacta como en el ejemplo dado. Creo que no estoy codificando las ecuaciones correctamente. Abajo está mi código

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

Mi salida para la desviación estándar es

0.1972484 0.1907743 0.1790021 0.1789577 0.1789325 0.2039248 0.1791031 0.2405984 0.1849784 0.1793203

y las existencias en el momento 1 son

50.36042 53.11321 64.32868 64.66536 66.02844 48.05690 63.80079 39.37478 56.44494 63.03220

mientras que en el ejemplo el autor tiene una salida diferente

Agradecería la ayuda.

Answer 1

Para ser honesto, no entiendo del PDF enlazado lo que lambda es o debe hacer, pero me las arreglé para igualar sus precios de las acciones con sólo dos cambios menores:

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

El tipo de interés y la volatilidad son anualizados, pero las simulaciones se basan en días. A menudo, la anualización se hace en base a 250 o 252 días, en este caso se utilizan 365 días. El código pasa a ser

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