<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/problems-in-estimating-garch-parameters-in-r/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://ntguardian.wordpress.com/2017/11/02/problems-estimating-garch-parameters-r/"> R – Curtis Miller's Personal Website</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p>These days my research focuses on change point detection methods. These are statistical tests and procedures to detect a structural change in a sequence of data. An early example, from quality control, is detecting whether a machine became uncalibrated when producing a widget. There may be some measurement of interest, such as the diameter of a ball bearing, that we observe. The machine produces these widgets in sequence. Under the null hypothesis, the ball bearing’s mean diameter does not change, while under the alternative, at some unkown point in the manufacturing process the machine became uncalibrated and the mean diameter of the ball bearings changed. The test then decides between these two hypotheses.</p>
<p><span></span></p>
<p>These types of test matter to economists and financial sector workers as well, particularly for forecasting. Once again we have a sequence of data indexed by time; my preferred example is price of a stock, which people can instantly recognize as a time series given how common time series graphs for stocks are, but there are many more datasets such as a state’s GDP or the unemployment rate. Economists want to forecast these quantities using past data and statistics. One of the assumptions the statistical methods makes is that the series being forecasted is stationary: the data was generated by one process with a single mean, autocorrelation, distribution, etc. This assumption isn’t always tested yet it is critical to successful forecasting. Tests for structural change check this assumption, and if it turns out to be false, the forecaster may need to divide up their dataset when training their models.</p>
<p>I have written about these tests <a rel="nofollow" href="https://ntguardian.wordpress.com/2016/09/12/wages-detach-productivity-1973/" target="_blank">before</a>, introducing the CUSUM statistic, one of the most popular statistics for detecting structural change. My advisor and a former Ph.D. student of his (currently a professor at the University of Waterloo, Greg Rice) developed a new test statistic that better detects structural changes that occur early or late in the dataset (imagine the machine producing widgets became uncalibrated just barely, and only the last dozen of the hundred widgets in the sample were affected). We’re in the process of making revisions requested by a journal to whom we submitted our paper, one of the revisions being a better example application (we initially worked with the wage/productivity data I discussed in the aforementioned blog post; the reviewers complained that these variables are codetermined so its nonsense to regress one on the other, a complaint I disagree with but I won’t plant my flag on to defend).</p>
<p>We were hoping to apply a version of our test to detecting structural change in GARCH models, a common model in financial time series. To my knowledge the “state of the art” R package for GARCH model estimation and inference (along with other work) is <a rel="nofollow" href="https://cran.r-project.org/web/packages/fGarch/index.html" target="_blank"><strong>fGarch</strong></a>; in particular, the function <code>garchFit()</code> is used for estimating GARCH models from data. When we tried to use this function in our test, though, we were given obviously bad numbers (we had already done simulation studies to know what behavior to expect). The null hypothesis of no change was soundly rejected on simulated sequences where it was true. I never saw the test fail to reject the null hypothesis, even though the null hypothesis was always true. This was the case even for sample sizes of 10,000, hardly a small sample.</p>
<p>We thought the problem might lie with the estimation of the covariance matrix of the parameter estimates, and I painstakingly derived and programmed functions to get this matrix not using numerical differentiation procedures, yet this did not stop the bad behavior. Eventually my advisor and I last Wednesday played with <code>garchFit()</code> and decided that the function is to blame. The behavior of the function on simulated data is so erratic when <em>estimating parameters</em> (not necessarily the covariance matrix as we initially thought, though it’s likely polluted as well) the function is basically useless, to my knowledge.</p>
<p>This function should be well-known and it’s certainly possible that the problem lies with me, not <strong>fGarch</strong> (or perhaps there’s better packages out there). This strikes me as a function of such importance I should share my findings. In this article I show a series of numerical experiments demonstrating <code>garchFit()</code>‘s pathological behavior.</p>
<h2>Basics on GARCH Models</h2>
<p>The <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> model is a time series model often used to model the volatility of financial instrument returns, such as the returns from stocks. Let <img alt="epsilon_t" title="epsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> represent the <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> process. This could represent the deviations in the returns of, say, a stock. The <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> model (without a mean parameter) is defined recursively as:</p>
<p><img alt="epsilon_t = sigma_t eta_t" title="epsilon_t = sigma_t eta_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t+%3D+%5Csigma_t+%5Ceta_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="sigma_t^2 = omega + alpha epsilon_{t - 1}^2 + beta sigma_{t - 1}^2" title="sigma_t^2 = omega + alpha epsilon_{t - 1}^2 + beta sigma_{t - 1}^2" src="https://s0.wp.com/latex.php?latex=%5Csigma_t%5E2+%3D+%5Comega+%2B+%5Calpha+%5Cepsilon_%7Bt+-+1%7D%5E2+%2B+%5Cbeta+%5Csigma_%7Bt+-+1%7D%5E2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="sigma_t" title="sigma_t" src="https://s0.wp.com/latex.php?latex=%5Csigma_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the conditional standard deviation of the process, also known as the conditional volatility, and <img alt="eta_t" title="eta_t" src="https://s0.wp.com/latex.php?latex=%5Ceta_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is a random process.</p>
<p>People who follow finance<sup><a rel="nofollow" href="https://ntguardian.wordpress.com/2017/11/02/problems-estimating-garch-parameters-r/#fn-3275-1" target="_blank">1</a></sup> noticed that returns to financial instruments (such as stocks or mutual funds) exhibit behavior known as volatility clustering. Some periods a financial instrument is relatively docile; there are not dramatic market movements. In others an instrument’s price can fluctuate greatly, and these periods are not one-off single-day movements but can last for a period of time. GARCH models were developed to model volatility clustering.</p>
<p>It is believed by some that even if a stock’s daily movement is essentially unforecastable (a stock is equally likely to over- or under-perform on any given day), the volatility <em>is</em> forecastable. Even for those who don’t have the hubris to believe anything about future returns can be forecasted these models are important. For example if one uses the model <img alt="_t = alpha + beta R_{M,t} + epsilon_t" title="_t = alpha + beta R_{M,t} + epsilon_t" src="https://s0.wp.com/latex.php?latex=_t+%3D+%5Calpha+%2B+%5Cbeta+R_%7BM%2Ct%7D+%2B+%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> to estimate the beta statistic for a stock (where <img alt="_t" title="_t" src="https://s0.wp.com/latex.php?latex=_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the stock’s return at time $latex $ <img alt="_{M,t}" title="_{M,t}" src="https://s0.wp.com/latex.php?latex=_%7BM%2Ct%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the market return, and <img alt="epsilon_t" title="epsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is “random noise”), there is a good chance that <img alt="epsilon_t" title="epsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is not an <em>i.i.d</em> sequence of random numbers (as is commonly assumed in other statistical contexts) but actually a GARCH sequence. The modeller would then want to know the behavior of her estimates in such a situation. Thus GARCH models are considered important. In fact, the volatility clustering behavior I just described is sometimes described as “GARCH behavior”, since it appears frequently and GARCH models are a frequent tool of choice to address them. (The acronym GARCH stands for generalized autoregressive conditional heteroskedasticity, which is statistics-speak for changing, time-dependent volatility.)</p>
<p><img alt="eta_t" title="eta_t" src="https://s0.wp.com/latex.php?latex=%5Ceta_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> can be any random process but a frequent choice is to use a sequence of <em>i.i.d</em> standard Normal random variables. Here <img alt="eta_t" title="eta_t" src="https://s0.wp.com/latex.php?latex=%5Ceta_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is the only source of randomness in the model. In order for a <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> process to have a stationary solution, we must require that <img alt="alpha + beta  0" title="alpha + beta  0" src="https://s0.wp.com/latex.php?latex=%5Calpha+%2B+%5Cbeta++0&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">). In this case the process has a long-run variance of <img alt="frac{omega}{1 - alpha - beta}" title="frac{omega}{1 - alpha - beta}" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B%5Comega%7D%7B1+-+%5Calpha+-+%5Cbeta%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">.</p>
<h2>Estimating GARCH Parameters</h2>
<p>The process I wrote down above is an <em>infinite</em> process; the index $latex $ can extend to negative numbers and beyond. Obviously in practice we don’t observe infinite sequences so if we want to work with <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> models in practice we need to consider a similar sequence:</p>
<p><img alt="epsilon_t = tilde{sigma}_t eta_t" title="epsilon_t = tilde{sigma}_t eta_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t+%3D+%5Ctilde%7B%5Csigma%7D_t+%5Ceta_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p><img alt="tilde{sigma}_t^2 = omega + alpha epsilon_{t - 1}^2 + beta tilde{sigma}_{t - 1}^2" title="tilde{sigma}_t^2 = omega + alpha epsilon_{t - 1}^2 + beta tilde{sigma}_{t - 1}^2" src="https://s0.wp.com/latex.php?latex=%5Ctilde%7B%5Csigma%7D_t%5E2+%3D+%5Comega+%2B+%5Calpha+%5Cepsilon_%7Bt+-+1%7D%5E2+%2B+%5Cbeta+%5Ctilde%7B%5Csigma%7D_%7Bt+-+1%7D%5E2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>Below is the new sequence’s secret sauce:</p>
<p><img alt="tilde{sigma}_0^2 = epsilon_0^2 = frac{omega}{1 - alpha - beta}" title="tilde{sigma}_0^2 = epsilon_0^2 = frac{omega}{1 - alpha - beta}" src="https://s0.wp.com/latex.php?latex=%5Ctilde%7B%5Csigma%7D_0%5E2+%3D+%5Cepsilon_0%5E2+%3D+%5Cfrac%7B%5Comega%7D%7B1+-+%5Calpha+-+%5Cbeta%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>We <em>choose</em> an initial value for this sequence (the theoretical <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> sequence described earlier does not have an initial value)! This sequence strongly resembles the theoretical sequence but it is observable in its entirity, and it can be shown that parameters estimated using this sequence closely approximate those of the theoretical, infinite <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> process.</p>
<p>Naturally one of the most important tasks for these processes is estimating their parameters; for the <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> process, these are <img alt="omega" title="omega" src="https://s0.wp.com/latex.php?latex=%5Comega&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, <img alt="alpha" title="alpha" src="https://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, and <img alt="beta" title="beta" src="https://s0.wp.com/latex.php?latex=%5Cbeta&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. A basic approach is to find the quasi-maximum likelihood estimation (QMLE) estimates. Let’s assume that we have $latex $ observations from our process. In QMLE, we work with the condisional distribution of <img alt="epsilon_t" title="epsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> when assuming <img alt="eta_t" title="eta_t" src="https://s0.wp.com/latex.php?latex=%5Ceta_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> follows a standard normal distribution (that is, <img alt="eta_t sim N(0, 1)" title="eta_t sim N(0, 1)" src="https://s0.wp.com/latex.php?latex=%5Ceta_t+%5Csim+N%280%2C+1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">). We assume that the entire history of the process up to time $latex $ is known; this implies that <img alt="tilde{sigma}_t" title="tilde{sigma}_t" src="https://s0.wp.com/latex.php?latex=%5Ctilde%7B%5Csigma%7D_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> is known as well (in fact all we needed to know was the values of the process at time <img alt="- 1" title="- 1" src="https://s0.wp.com/latex.php?latex=-+1&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, but I digress). In that case we have <img alt="epsilon_t sim N(0, tilde{sigma}_t^2)" title="epsilon_t sim N(0, tilde{sigma}_t^2)" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t+%5Csim+N%280%2C+%5Ctilde%7B%5Csigma%7D_t%5E2%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. Let <img alt="(x | mathscr{f}_t) = f(x | tilde{sigma}_t^2)" title="(x | mathscr{f}_t) = f(x | tilde{sigma}_t^2)" src="https://s0.wp.com/latex.php?latex=%28x+%7C+%5Cmathscr%7Bf%7D_t%29+%3D+f%28x+%7C+%5Ctilde%7B%5Csigma%7D_t%5E2%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> be the conditional distribution of <img alt="epsilon_t" title="epsilon_t" src="https://s0.wp.com/latex.php?latex=%5Cepsilon_t&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> (so <img alt="(t | tilde{sigma}_t^2) = frac{1}{sqrt{2 pi tilde{sigma}_t^2}} expleft(-frac{x^2}{2 tilde{sigma}_t^2}right)" title="(t | tilde{sigma}_t^2) = frac{1}{sqrt{2 pi tilde{sigma}_t^2}} expleft(-frac{x^2}{2 tilde{sigma}_t^2}right)" src="https://s0.wp.com/latex.php?latex=%28t+%7C+%5Ctilde%7B%5Csigma%7D_t%5E2%29+%3D+%5Cfrac%7B1%7D%7B%5Csqrt%7B2+%5Cpi+%5Ctilde%7B%5Csigma%7D_t%5E2%7D%7D+%5Cexp%5Cleft%28-%5Cfrac%7Bx%5E2%7D%7B2+%5Ctilde%7B%5Csigma%7D_t%5E2%7D%5Cright%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">). The quasi-likelihood equation is then</p>
<p><img alt="mathscr{L}_T(epsilon_1, ..., epsilon_T) = prod_{t = 1}^T f(epsilon_t | tilde{sigma}_t^2) = (2 pi)^{-n/2} prod_{t = 1}^T tilde{sigma}_t^{-1} expleft(-frac{1}{2}sum_{t = 1}^T frac{epsilon_t^2}{tilde{sigma}_t^2} right)" title="mathscr{L}_T(epsilon_1, ..., epsilon_T) = prod_{t = 1}^T f(epsilon_t | tilde{sigma}_t^2) = (2 pi)^{-n/2} prod_{t = 1}^T tilde{sigma}_t^{-1} expleft(-frac{1}{2}sum_{t = 1}^T frac{epsilon_t^2}{tilde{sigma}_t^2} right)" src="https://s0.wp.com/latex.php?latex=%5Cmathscr%7BL%7D_T%28%5Cepsilon_1%2C+...%2C+%5Cepsilon_T%29+%3D+%5Cprod_%7Bt+%3D+1%7D%5ET+f%28%5Cepsilon_t+%7C+%5Ctilde%7B%5Csigma%7D_t%5E2%29+%3D+%282+%5Cpi%29%5E%7B-n%2F2%7D+%5Cprod_%7Bt+%3D+1%7D%5ET+%5Ctilde%7B%5Csigma%7D_t%5E%7B-1%7D+%5Cexp%5Cleft%28-%5Cfrac%7B1%7D%7B2%7D%5Csum_%7Bt+%3D+1%7D%5ET+%5Cfrac%7B%5Cepsilon_t%5E2%7D%7B%5Ctilde%7B%5Csigma%7D_t%5E2%7D+%5Cright%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>Like most likelihood methods, rather than optimize the quasi-likelihood function directly, statisticians try to optimize the log-likelihood, <img alt="logleft(mathscr{L}_T(epsilon_1, ..., epsilon_T)right)" title="logleft(mathscr{L}_T(epsilon_1, ..., epsilon_T)right)" src="https://s0.wp.com/latex.php?latex=%5Clog%5Cleft%28%5Cmathscr%7BL%7D_T%28%5Cepsilon_1%2C+...%2C+%5Cepsilon_T%29%5Cright%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, and after some work it’s not hard to see this is equivalent to minimizing</p>
<p><img alt="sum_{t = 1}^T left( log(tilde{sigma}_t^2) + frac{epsilon_t^2}{tilde{sigma}_t^2} right)" title="sum_{t = 1}^T left( log(tilde{sigma}_t^2) + frac{epsilon_t^2}{tilde{sigma}_t^2} right)" src="https://s0.wp.com/latex.php?latex=%5Csum_%7Bt+%3D+1%7D%5ET+%5Cleft%28+%5Clog%28%5Ctilde%7B%5Csigma%7D_t%5E2%29+%2B+%5Cfrac%7B%5Cepsilon_t%5E2%7D%7B%5Ctilde%7B%5Csigma%7D_t%5E2%7D+%5Cright%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"></p>
<p>Note that <img alt="omega" title="omega" src="https://s0.wp.com/latex.php?latex=%5Comega&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, <img alt="alpha" title="alpha" src="https://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, and <img alt="beta" title="beta" src="https://s0.wp.com/latex.php?latex=%5Cbeta&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> are involved in this quantity through <img alt="tilde{sigma}_t^2" title="tilde{sigma}_t^2" src="https://s0.wp.com/latex.php?latex=%5Ctilde%7B%5Csigma%7D_t%5E2&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">. There is no closed form solution for the parameters that minimize this quantity. This means that numerical optimization techniques must be applied to find the parameters.</p>
<p>It can be shown that the estimators for the parameters <img alt="omega" title="omega" src="https://s0.wp.com/latex.php?latex=%5Comega&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, <img alt="alpha" title="alpha" src="https://s0.wp.com/latex.php?latex=%5Calpha&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, and <img alt="beta" title="beta" src="https://s0.wp.com/latex.php?latex=%5Cbeta&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex">, when computed this way, are consistent (meaning that asymptotically they approach their true values, in the sense that they converge in probability) and follow a Gaussian distribution asymptotically.<sup><a rel="nofollow" href="https://ntguardian.wordpress.com/2017/11/02/problems-estimating-garch-parameters-r/#fn-3275-2" target="_blank">2</a></sup> These are properties that we associate with the sample mean, and while we might be optimistic that the rate of convergence of these estimators is as good as the rate of convergence of the sample mean, we may expect comparable asymptotic behavior.</p>
<p>Ideally, the parameters should behave like the process illustrated below.</p>
<pre>
library(ggplot2)

x &lt;- rnorm(1000, sd = 1/3)
df &lt;- t(sapply(50:1000, function(t) {
  return(c(&quot;mean&quot; = mean(x[1:t]), &quot;mean.se&quot; = sd(x[1:t])/sqrt(t)))
}))
df &lt;- as.data.frame(df)
df$t &lt;- 50:1000
ggplot(df, aes(x = t, y = mean)) +
  geom_line() +
  geom_ribbon(aes(x = t, ymin = mean - 2 * mean.se, ymax = mean + 2 * mean.se),
              color = &quot;grey&quot;, alpha = 0.5) +
  geom_hline(color = &quot;blue&quot;, yintercept = 0) +
  coord_cartesian(ylim = c(-0.5, 0.5))
</pre>
<p><img alt="" src="https://ntguardian.files.wordpress.com/2017/11/unnamed-chunk-2-1.png?w=456" data-recalc-dims="1"></p>
<h2>Behavior of Estimates by <strong>fGarch</strong></h2>
<p>Before continuing let’s generate a <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> sequence. Throughout this article I work with processes where all parameters are equal to 0.2. Notice that for a <img alt="text{GARCH}(1,1)" title="text{GARCH}(1,1)" src="https://s0.wp.com/latex.php?latex=%5Ctext%7BGARCH%7D%281%2C1%29&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> process the long-run variance will be <img alt="frac{0.2}{1 - 0.2 - 0.2} = frac{1}{3}" title="frac{0.2}{1 - 0.2 - 0.2} = frac{1}{3}" src="https://s0.wp.com/latex.php?latex=%5Cfrac%7B0.2%7D%7B1+-+0.2+-+0.2%7D+%3D+%5Cfrac%7B1%7D%7B3%7D&bg=ffffff&%23038;fg=444444&%23038;s=0" class="latex"> with this choice.</p>
<pre>
set.seed(110117)
library(fGarch)
x &lt;- garchSim(garchSpec(model = list(&quot;alpha&quot; = 0.2, &quot;beta&quot; = 0.2, &quot;omega&quot; = 0.2)),
              n.start = 1000, n = 1000)
plot(x)
</pre>
<p><img alt="" src="https://ntguardian.files.wordpress.com/2017/11/unnamed-chunk-3-1.png?w=456" data-recalc-dims="1"></p>
<p>Let’s see the parameters that the <strong>fGarch</strong> function <code>garchFit()</code> uses.</p>
<pre>
args(garchFit)
</pre>
<pre>
## function (formula = ~garch(1, 1), data = dem2gbp, init.rec = c(&quot;mci&quot;,
##     &quot;uev&quot;), delta = 2, skew = 1, shape = 4, cond.dist = c(&quot;norm&quot;,
##     &quot;snorm&quot;, &quot;ged&quot;, &quot;sged&quot;, &quot;std&quot;, &quot;sstd&quot;, &quot;snig&quot;, &quot;QMLE&quot;), include.mean = TRUE,
##     include.delta = NULL, include.skew = NULL, include.shape = NULL,
##     leverage = NULL, trace = TRUE, algorithm = c(&quot;nlminb&quot;, &quot;lbfgsb&quot;,
##         &quot;nlminb+nm&quot;, &quot;lbfgsb+nm&quot;), hessian = c(&quot;ropt&quot;, &quot;rcd&quot;),
##     control = list(), title = NULL, description = NULL, ...)
## NULL
</pre>
<p>The function provides a few options for distribution to maximize (<code>cond.dist</code>) and algorithm to use for optimization (<code>algorithm</code>). Here I will always choose <code>cond.dist = QMLE</code>, unless otherwise stated, to instruct the function to use QMLE estimators.</p>
<p>Here’s a single pass.</p>
<pre>
garchFit(data = x, cond.dist = &quot;QMLE&quot;, include.mean = FALSE)
</pre>
<pre>
##
## Series Initialization:
##  ARMA Model:                arma
##  Formula Mean:              ~ arma(0, 0)
##  GARCH Model:               garch
##  Formula Variance:          ~ garch(1, 1)
##  ARMA Order:                0 0
##  Max ARMA Order:            0
##  GARCH Order:               1 1
##  Max GARCH Order:           1
##  Maximum Order:             1
##  Conditional Dist:          QMLE
##  h.start:                   2
##  llh.start:                 1
##  Length of Series:          1000
##  Recursion Init:            mci
##  Series Scale:              0.5320977
##
## Parameter Initialization:
##  Initial Parameters:          $params
##  Limits of Transformations:   $U, $V
##  Which Parameters are Fixed?  $includes
##  Parameter Matrix:
##                      U          V params includes
##     mu     -0.15640604   0.156406    0.0    FALSE
##     omega   0.00000100 100.000000    0.1     TRUE
##     alpha1  0.00000001   1.000000    0.1     TRUE
##     gamma1 -0.99999999   1.000000    0.1    FALSE
##     beta1   0.00000001   1.000000    0.8     TRUE
##     delta   0.00000000   2.000000    2.0    FALSE
##     skew    0.10000000  10.000000    1.0    FALSE
##     shape   1.00000000  10.000000    4.0    FALSE
##  Index List of Parameters to be Optimized:
##  omega alpha1  beta1
##      2      3      5
##  Persistence:                  0.9
##
##
## --- START OF TRACE ---
## Selected Algorithm: nlminb
##
## R coded nlminb Solver:
##
##   0:     1419.0152: 0.100000 0.100000 0.800000
##   1:     1418.6616: 0.108486 0.0998447 0.804683
##   2:     1417.7139: 0.109746 0.0909961 0.800931
##   3:     1416.7807: 0.124977 0.0795152 0.804400
##   4:     1416.7215: 0.141355 0.0446605 0.799891
##   5:     1415.5139: 0.158059 0.0527601 0.794304
##   6:     1415.2330: 0.166344 0.0561552 0.777108
##   7:     1415.0415: 0.195230 0.0637737 0.743465
##   8:     1415.0031: 0.200862 0.0576220 0.740088
##   9:     1414.9585: 0.205990 0.0671331 0.724721
##  10:     1414.9298: 0.219985 0.0713468 0.712919
##  11:     1414.8226: 0.230628 0.0728325 0.697511
##  12:     1414.4689: 0.325750 0.0940514 0.583114
##  13:     1413.4560: 0.581449 0.143094 0.281070
##  14:     1413.2804: 0.659173 0.157127 0.189282
##  15:     1413.2136: 0.697840 0.155964 0.150319
##  16:     1413.1467: 0.720870 0.142550 0.137645
##  17:     1413.1416: 0.726527 0.138146 0.135966
##  18:     1413.1407: 0.728384 0.137960 0.134768
##  19:     1413.1392: 0.731725 0.138321 0.132991
##  20:     1413.1392: 0.731146 0.138558 0.133590
##  21:     1413.1392: 0.730849 0.138621 0.133850
##  22:     1413.1392: 0.730826 0.138622 0.133869
##
## Final Estimate of the Negative LLH:
##  LLH:  782.211    norm LLH:  0.782211
##     omega    alpha1     beta1
## 0.2069173 0.1386221 0.1338686
##
## R-optimhess Difference Approximated Hessian Matrix:
##            omega     alpha1      beta1
## omega  -8858.897 -1839.6144 -2491.9827
## alpha1 -1839.614  -782.8005  -531.7393
## beta1  -2491.983  -531.7393  -729.7246
## attr(,&quot;time&quot;)
## Time difference of 0.04132652 secs
##
## --- END OF TRACE ---
##
##
## Time to Estimate Parameters:
##  Time difference of 0.3866439 secs
</pre>
<pre>
##
## Title:
##  GARCH Modelling
##
## Call:
##  garchFit(data = x, cond.dist = &quot;QMLE&quot;, include.mean = FALSE)
##
## Mean and Variance Equation:
##  data ~ garch(1, 1)
## </pre>