GARCH(1,1)-M MLE optimization with fmincon in R

External question — Quantitative Finance Stack Exchange Author: SlavicDoomer Original post: https://quant.stackexchange.com/questions/51645 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. I've searched thru dozens of papers and did not find in any of them satisfying and enough theoretical answers to my concerns. So I've combined everything what I found below. Please indicate if my understanding of the topic is proper and correct me if necessary. Beside theoretical side I've found also an obstacle in practical implementation. Duan (1995) in his paper developed model for European options pricing with GARCH model. He introduced locally risk neutral measure (LRNVR) $\mathbb{Q} $ which is equivalent to physical market measure $\mathbb{P} $ . Since price for either call or put option has no analytical solution in his framework it is necessary to run Monte Carlo simulations. My first concern for which I did not found explicit answer is the following: (1) We have a time serie with some market data and we fit parameters of the model to this time serie under the physical measure $ \mathbb{P} $ . (2) Then we use estimated parameters to the transformed process under the LRNVR measure $\mathbb{Q} $ and we run Monte Carlo simulations to estimate price of an option. Since MC simulation under measure $\mathbb{Q} $ is not a problem for me, let's now focus on theoretical derivation of MLE under physical measure $\mathbb{P} $. Suppose we have a sample of $T$ log-returns of some financial asset. Let $X = (X_1, \dots , X_T) $ denote our sample and assume that $t$ -th log-return follows GARCH(1,1)-M process proposed by Duan (1995): $$ X_t = \ln \frac{S_t}{S_{t-1}} = r + \lambda \sqrt{h_t} - 0.5h_t + \sqrt{h_t} z_t $$ under physical measure $\mathbb{P} $ , where $z_t \overset{iid}{=} \mathcal{N} (0,1) $ and $h_t = \omega + \alpha h_{t-1} z_{t-1}^2 + \beta h_{t-1} $ . Here we assume that $ \omega >0 $ and $ \alpha, \beta \geq 0 $ . We would also require $ \alpha + \beta < 1$ to ensure stationarity. Parameter $r$ is a market risk-free interest rate (known parameter) and $\lambda$ is a risk-premium associated with a given financial asset (parameter to be estimated). Since we have $ X_t \sim \mathcal{N} (r + \lambda \sqrt{h_t} - 0.5h_t, h_t )$ , then likelihood function for $t$ -th observation is: $$ l_t (X_t ; \theta ) = \frac{1}{\sqrt{2 \pi h_t}} \exp \left( - \frac{ ( X_t - r - \lambda \sqrt{h_t} + \frac{1}{2} h_t )^2 }{2h_t} \right) \text{,} $$ where $\theta = (\omega, \alpha, \beta, \lambda ) $ is a vector of parameters to be estimated. Likelihood function for a vector $X$ is: $$ \mathcal{L} (X; \theta ) = \prod_{t=1}^T l_t(X_t; \theta) \text{.} $$ Since it is easier to compute natural logs, then: $$ \ln \mathcal{L} (X; \theta ) = -\frac{T}{2} \ln \left( 2 \pi \right) - \frac{1}{2} \sum_{t=1}^T \ln \left( h_t \right) - \frac{1}{2} \sum_{t=1}^T \frac{ ( X_t - r - \lambda \sqrt{h_t} + \frac{1}{2} h_t )^2 }{2h_t} $$ to be maximized. We are seeking for $ \hat{\theta} $ which is: $$ \hat{\theta} = \arg \max_{\theta} \ln \mathcal{L} (X; \theta ) $$ with constrains for $\omega , \alpha, \beta $ as indicated above. Now let's move to the practical implementation of the above: Let our market data will be AAPL daily log-returns for the period 2016-2019 (or any other data, because the following problem does not disappear with the change of source data). Suppose risk-free interest rate is $r=0$ . As initial variance $h_1$ we assume variance of our sample, i.e. $h_1 = Var(X) $ . Function to be minimized is $ - \ln \mathcal{L} (X; \theta ) $ and is defined as follows: loglike <- function(params, log_returns){ omega <- params[1] alpha <- params[2] beta <- params[3] lambda <- params[4] bigT <- length(log_returns) h <- c(var(log_returns)) for (i in 2:bigT) { h[i] <- omega + alpha * (log_returns[i-1] - lambda * sqrt(h[i-1]) + 0.5 * h[i-1] )^2 + beta * h[i-1] } likelihood <- 0.5 * bigT * log(2*pi) + 0.5 * sum(log(h) + ((log_returns - lambda * sqrt(h) + 0.5 * h)^2) / h ) return(likelihood) } Starting parameters and constrains for optimization problem are as follows: params <- rep(0.01, 4) lb <- c(0,0,0, -Inf) A = cbind(0, 1, 1, 0) b = 1 And if I try to use fmincon from pracma package for optimization I have the following error and warnings: > fmincon(loglike, log_returns = log_returns, x0 = params, lb = lb, A = A, b = b) Error in if (f < 0) { : missing value where TRUE/FALSE needed In addition: Warning messages: 1: In sqrt(h[i - 1]) : NaNs produced 2: In log(h) : NaNs produced 3: In sqrt(h) : NaNs produced First of all I do not know why I get these warnings. According to the objective function defined above $ h \geq 0 $ so taking square roots and natural logs should not produce NaNs. Secondly I do not understand the error returned by fmincon . What is wrong with my objective function?

External answer — Quantitative Finance Stack Exchange Author: Andreas Original post: https://quant.stackexchange.com/a/53269 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. There are two things i recommend you to try, maybe that already fixes your problem. First, drop the constant term $-\frac{T}{2}ln\left ( 2\pi \right )$ to eliminate one possible source for your error. (this part it is not needed for the optimization anyways) Second, sometimes errors arise due to numerical issues when optimization routines do not strictly stick to the limits set by the user (even if they should). You can try setting your limit to an arbitrarily small positive number.

External answer — Quantitative Finance Stack Exchange Author: QuantCalc.net Original post: https://quant.stackexchange.com/a/85298 License: CC BY-SA 4.0 — https://creativecommons.org/licenses/by-sa/4.0/ Adaptation: HTML converted to plain text; contact email addresses removed. The main issue is almost certainly within the loop that calculates the conditional variance, $h_i$ . If $h_i=0$ , you have division by zero issue and $\ln(0)$ issue in the likelihood function. The simplest way to enforce this is to constrain the parameters: ω>0 (Omega) α≥0 (Alpha) β≥0 (Beta) You must ensure ω is strictly positive by using a small positive number (e.g., 1e−6) in the lb (lower bound) argument. lb <- c(1e-6, 0, 0, -Inf) ​