Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_tsa_uni_garch_exp_estim (g13fg)

## Purpose

nag_tsa_uni_garch_exp_estim (g13fg) estimates the parameters of a univariate regression-exponential GARCH(p,q)$\text{GARCH}\left(p,q\right)$ process (see Engle and Ng (1993)).

## Syntax

[theta, se, sc, covr, hp, et, ht, lgf, ifail] = g13fg(dist, yt, x, ip, iq, nreg, mn, theta, hp, copts, maxit, tol, 'num', num, 'npar', npar)
[theta, se, sc, covr, hp, et, ht, lgf, ifail] = nag_tsa_uni_garch_exp_estim(dist, yt, x, ip, iq, nreg, mn, theta, hp, copts, maxit, tol, 'num', num, 'npar', npar)

## Description

A univariate regression-exponential GARCH(p,q)$\text{GARCH}\left(p,q\right)$ process, with q$q$ coefficients αi${\alpha }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$, q$q$ coefficients φi${\varphi }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$, p$p$ coefficients, βi${\beta }_{\mathit{i}}$, for i = 1,2,,p$\mathit{i}=1,2,\dots ,p$, and k$k$ linear regression coefficients bi${b}_{\mathit{i}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$, can be represented by:
yt = bo + xtT b + εt
 q q p ln(ht) = α0 + ∑ αizt − i + ∑ φi (|zt − i| − E[|zt − i|]) + ∑ βiln(ht − i),  t = 1,2, … ,T i = 1 i = 1 i = 1
$yt = bo + xtT b + εt ln(ht) = α0+ ∑ i=1 q αi zt-i + ∑ i=1 q ϕi ( | zt-i | - E [ | zt-i | ] ) + ∑ i=1 p βi ln ( ht-i ) , t=1,2,…,T$
(1)
where zt = (εt)/(sqrt(ht)) ${z}_{t}=\frac{{\epsilon }_{t}}{\sqrt{{h}_{t}}}$, E[|zti|]$E\left[|{z}_{t-i}|\right]$ denotes the expected value of |zti|$|{z}_{t-i}|$ and εtψt1 = N(0,ht)${\epsilon }_{t}\mid {\psi }_{t-1}=N\left(0,{h}_{t}\right)$ or εtψt1 = St(df,ht)${\epsilon }_{t}\mid {\psi }_{t-1}={S}_{t}\left(\mathit{df},{h}_{t}\right)$. Here St${S}_{t}$ is a standardized Student's t$t$-distribution with df$\mathit{df}$ degrees of freedom and variance ht${h}_{t}$, T$T$ is the number of terms in the sequence, yt${y}_{t}$ denotes the endogenous variables, xt${x}_{t}$ the exogenous variables, bo${b}_{o}$ the regression mean, b$b$ the regression coefficients, εt${\epsilon }_{t}$ the residuals, ht${h}_{t}$ the conditional variance, df$\mathit{df}$ the number of degrees of freedom of the Student's t$t$-distribution, and ψt${\psi }_{t}$ the set of all information up to time t$t$.
nag_tsa_uni_garch_exp_estim (g13fg) provides an estimate θ̂$\stackrel{^}{\theta }$, for the vector θ = (bo,bT,ωT) $\theta =\left({b}_{o},{b}^{\mathrm{T}},{\omega }^{\mathrm{T}}\right)$ where bT = (b1,,bk) ${b}^{\mathrm{T}}=\left({b}_{1},\dots ,{b}_{k}\right)$, ωT = (α0,α1,,αq,φ1,,φq,β1,,βp,γ) ${\omega }^{\mathrm{T}}=\left({\alpha }_{0},{\alpha }_{1},\dots ,{\alpha }_{q},{\varphi }_{1},\dots ,{\varphi }_{q},{\beta }_{1},\dots ,{\beta }_{p},\gamma \right)$ when dist = 'N'${\mathbf{dist}}=\text{'N'}$, and ωT = (α0,α1,,αq,φ1,,φq,β1,,βp,γ,df) ${\omega }^{\mathrm{T}}=\left({\alpha }_{0},{\alpha }_{1},\dots ,{\alpha }_{q},{\varphi }_{1},\dots ,{\varphi }_{q},{\beta }_{1},\dots ,{\beta }_{p},\gamma ,\mathit{df}\right)$ when dist = 'T'${\mathbf{dist}}=\text{'T'}$.
mn, nreg can be used to simplify the GARCH(p,q)$\text{GARCH}\left(p,q\right)$ expression in (1) as follows:
No Regression and No Mean
• yt = εt${y}_{t}={\epsilon }_{t}$,
• mn = 0${\mathbf{mn}}=0$,
• nreg = 0${\mathbf{nreg}}=0$ and
• θ$\theta$ is a (2 × q + p + 1)$\left(2×q+p+1\right)$ vector when dist = 'N'${\mathbf{dist}}=\text{'N'}$, and a (2 × q + p + 2)$\left(2×q+p+2\right)$ vector, when dist = 'T'${\mathbf{dist}}=\text{'T'}$.
No Regression
• yt = bo + εt${y}_{t}={b}_{o}+{\epsilon }_{t}$,
• mn = 1${\mathbf{mn}}=1$,
• nreg = 0${\mathbf{nreg}}=0$ and
• θ$\theta$ is a (2 × q + p + 2)$\left(2×q+p+2\right)$ vector when dist = 'N'${\mathbf{dist}}=\text{'N'}$ and a (2 × q + p + 3)$\left(2×q+p+3\right)$ vector, when dist = 'T'${\mathbf{dist}}=\text{'T'}$.
Note:  if the yt = μ + εt${y}_{t}=\mu +{\epsilon }_{t}$, where μ$\mu$ is known (not to be estimated by nag_tsa_uni_garch_exp_estim (g13fg)) then (1) can be written as ytμ = εt${y}_{t}^{\mu }={\epsilon }_{t}$, where ytμ = ytμ${y}_{t}^{\mu }={y}_{t}-\mu$. This corresponds to the case No Regression and No Mean, with yt${y}_{t}$ replaced by ytμ${y}_{t}-\mu$.
No Mean
• yt = xtT b + εt ${y}_{t}={x}_{t}^{\mathrm{T}}b+{\epsilon }_{t}$,
• mn = 0${\mathbf{mn}}=0$,
• nreg = k${\mathbf{nreg}}=k$ and
• θ$\theta$ is a (2 × q + p + 1 + k)$\left(2×q+p+1+k\right)$ vector when dist = 'N'${\mathbf{dist}}=\text{'N'}$ and a (2 × q + p + 2 + k)$\left(2×q+p+2+k\right)$ vector, when dist = 'T'${\mathbf{dist}}=\text{'T'}$.

## References

Bollerslev T (1986) Generalised autoregressive conditional heteroskedasticity Journal of Econometrics 31 307–327
Engle R (1982) Autoregressive conditional heteroskedasticity with estimates of the variance of United Kingdom inflation Econometrica 50 987–1008
Engle R and Ng V (1993) Measuring and testing the impact of news on volatility Journal of Finance 48 1749–1777
Glosten L, Jagannathan R and Runkle D (1993) Relationship between the expected value and the volatility of nominal excess return on stocks Journal of Finance 48 1779–1801
Hamilton J (1994) Time Series Analysis Princeton University Press

## Parameters

### Compulsory Input Parameters

1:     dist – string (length ≥ 1)
The type of distribution to use for et${e}_{t}$.
dist = 'N'${\mathbf{dist}}=\text{'N'}$
A Normal distribution is used.
dist = 'T'${\mathbf{dist}}=\text{'T'}$
A Student's t$t$-distribution is used.
Constraint: dist = 'N'${\mathbf{dist}}=\text{'N'}$ or 'T'$\text{'T'}$.
2:     yt(num) – double array
num, the dimension of the array, must satisfy the constraint
• nummax (ip,iq)${\mathbf{num}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$
• ${\mathbf{num}}\ge {\mathbf{nreg}}+{\mathbf{mn}}$
• .
The sequence of observations, yt${y}_{\mathit{t}}$, for t = 1,2,,T$\mathit{t}=1,2,\dots ,T$.
3:     x(ldx, : $:$) – double array
The first dimension of the array x must be at least num${\mathbf{num}}$
The second dimension of the array must be at least nreg${\mathbf{nreg}}$
Row t$\mathit{t}$ of x must contain the time dependent exogenous vector xt ${x}_{\mathit{t}}$, where xtT = (xt1,,xtk) ${x}_{\mathit{t}}^{\mathrm{T}}=\left({x}_{\mathit{t}}^{1},\dots ,{x}_{\mathit{t}}^{k}\right)$, for t = 1,2,,T$\mathit{t}=1,2,\dots ,T$.
4:     ip – int64int32nag_int scalar
The number of coefficients, βi${\beta }_{\mathit{i}}$, for i = 1,2,,p$\mathit{i}=1,2,\dots ,p$.
Constraint: ip0${\mathbf{ip}}\ge 0$ (see also npar).
5:     iq – int64int32nag_int scalar
The number of coefficients, αi${\alpha }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$.
Constraint: iq1${\mathbf{iq}}\ge 1$ (see also npar).
6:     nreg – int64int32nag_int scalar
k$k$, the number of regression coefficients.
Constraint: nreg0${\mathbf{nreg}}\ge 0$ (see also npar).
7:     mn – int64int32nag_int scalar
If mn = 1${\mathbf{mn}}=1$, the mean term b0${b}_{0}$ will be included in the model.
Constraint: mn = 0${\mathbf{mn}}=0$ or 1$1$.
8:     theta(npar) – double array
npar, the dimension of the array, must satisfy the constraint npar < 20${\mathbf{npar}}<20$.
The initial parameter estimates for the vector θ$\theta$.
The first element must contain the coefficient αo${\alpha }_{o}$ and the next iq elements must contain the autoregressive coefficients αi${\alpha }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$.
The next iq elements contain the coefficients φi${\varphi }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$.
The next ip elements must contain the moving average coefficients βi${\beta }_{\mathit{i}}$, for i = 1,2,,p$\mathit{i}=1,2,\dots ,p$.
If dist = 'T'${\mathbf{dist}}=\text{'T'}$, the next element must contain an estimate for df$\mathit{df}$, the number of degrees of freedom of the Student's t$t$-distribution.
If mn = 1${\mathbf{mn}}=1$, the next element must contain the mean term bo${b}_{o}$.
If copts = false${\mathbf{copts}}=\mathbf{false}$, the remaining nreg elements are taken as initial estimates of the linear regression coefficients bi${b}_{\mathit{i}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
9:     hp – double scalar
If copts = false${\mathbf{copts}}=\mathbf{false}$ then hp is the value to be used for the pre-observed conditional variance, otherwise hp is not referenced.
10:   copts – logical scalar
If copts = true${\mathbf{copts}}=\mathbf{true}$, the function provides initial parameter estimates of the regression terms, otherwise these are provided by you.
11:   maxit – int64int32nag_int scalar
The maximum number of iterations to be used by the optimization function when estimating the GARCH(p,q)$\text{GARCH}\left(p,q\right)$ parameters.
Constraint: maxit > 0${\mathbf{maxit}}>0$.
12:   tol – double scalar
The tolerance to be used by the optimization function when estimating the GARCH(p,q)$\text{GARCH}\left(p,q\right)$ parameters.

### Optional Input Parameters

1:     num – int64int32nag_int scalar
Default: The dimension of the array yt and the first dimension of the array x. (An error is raised if these dimensions are not equal.)
T$T$, the number of terms in the sequence.
Constraints:
• nummax (ip,iq)${\mathbf{num}}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$;
• ${\mathbf{num}}\ge {\mathbf{nreg}}+{\mathbf{mn}}$.
2:     npar – int64int32nag_int scalar
Default: The dimension of the array theta.
The number of parameters to be included in the model. npar = 1 + 2 × iq + ip + mn + nreg${\mathbf{npar}}=1+2×{\mathbf{iq}}+{\mathbf{ip}}+{\mathbf{mn}}+{\mathbf{nreg}}$ when dist = 'N'${\mathbf{dist}}=\text{'N'}$ and npar = 2 + 2 × iq + ip + mn + nreg${\mathbf{npar}}=2+2×{\mathbf{iq}}+{\mathbf{ip}}+{\mathbf{mn}}+{\mathbf{nreg}}$ when dist = 'T'${\mathbf{dist}}=\text{'T'}$.
Constraint: npar < 20${\mathbf{npar}}<20$.

### Input Parameters Omitted from the MATLAB Interface

ldx ldcovr work lwork

### Output Parameters

1:     theta(npar) – double array
The estimated values θ̂$\stackrel{^}{\theta }$ for the vector θ$\theta$.
The first element contains the coefficient αo${\alpha }_{o}$ and the next iq elements contain the coefficients αi${\alpha }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$.
The next iq elements contain the coefficients φi${\varphi }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$.
The next ip elements are the moving average coefficients βi${\beta }_{\mathit{i}}$, for i = 1,2,,p$\mathit{i}=1,2,\dots ,p$.
If dist = 'T'${\mathbf{dist}}=\text{'T'}$, the next element contains an estimate for df$\mathit{df}$ then the number of degrees of freedom of the Student's t$t$-distribution.
If mn = 1${\mathbf{mn}}=1$, the next element contains an estimate for the mean term bo${b}_{o}$.
The final nreg elements are the estimated linear regression coefficients bi${b}_{\mathit{i}}$, for i = 1,2,,k$\mathit{i}=1,2,\dots ,k$.
2:     se(npar) – double array
The standard errors for θ̂$\stackrel{^}{\theta }$.
The first element contains the standard error for αo${\alpha }_{o}$ and the next iq elements contain the standard errors for αi${\alpha }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$. The next iq elements contain the standard errors for φi${\varphi }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$. The next ip elements are the standard errors for βj${\beta }_{\mathit{j}}$, for j = 1,2,,p$\mathit{j}=1,2,\dots ,p$.
If dist = 'T'${\mathbf{dist}}=\text{'T'}$, the next element contains the standard error for df$\mathit{df}$, the number of degrees of freedom of the Student's t$t$-distribution.
If mn = 1${\mathbf{mn}}=1$, the next element contains the standard error for bo${b}_{o}$.
The final nreg elements are the standard errors for bj${b}_{\mathit{j}}$, for j = 1,2,,k$\mathit{j}=1,2,\dots ,k$.
3:     sc(npar) – double array
The scores for θ̂$\stackrel{^}{\theta }$.
The first element contains the scores for αo${\alpha }_{o}$, the next iq elements contain the scores for αi${\alpha }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$, the next iq elements contain the scores for φi${\varphi }_{\mathit{i}}$, for i = 1,2,,q$\mathit{i}=1,2,\dots ,q$, the next ip elements are the scores for βj${\beta }_{\mathit{j}}$, for j = 1,2,,p$\mathit{j}=1,2,\dots ,p$.
If dist = 'T'${\mathbf{dist}}=\text{'T'}$, the next element contains the scores for df$\mathit{df}$, the number of degrees of freedom of the Student's t$t$-distribution.
If mn = 1${\mathbf{mn}}=1$, the next element contains the score for bo${b}_{o}$.
The final nreg elements are the scores for bj${b}_{\mathit{j}}$, for j = 1,2,,k$\mathit{j}=1,2,\dots ,k$.
4:     covr(ldcovr,npar) – double array
ldcovrnpar$\mathit{ldcovr}\ge {\mathbf{npar}}$.
The covariance matrix of the parameter estimates θ̂$\stackrel{^}{\theta }$, that is the inverse of the Fisher Information Matrix.
5:     hp – double scalar
If copts = true${\mathbf{copts}}=\mathbf{true}$ then hp is the estimated value of the pre-observed conditional variance.
6:     et(num) – double array
The estimated residuals, εt${\epsilon }_{\mathit{t}}$, for t = 1,2,,T$\mathit{t}=1,2,\dots ,T$.
7:     ht(num) – double array
The estimated conditional variances, ht${h}_{\mathit{t}}$, for t = 1,2,,T$\mathit{t}=1,2,\dots ,T$.
8:     lgf – double scalar
The value of the log-likelihood function at θ̂$\stackrel{^}{\theta }$.
9:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Note: nag_tsa_uni_garch_exp_estim (g13fg) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:
ifail = 1${\mathbf{ifail}}=1$
 On entry, nreg < 0${\mathbf{nreg}}<0$, or mn > 1${\mathbf{mn}}>1$, or mn < 0${\mathbf{mn}}<0$, or iq < 1${\mathbf{iq}}<1$, or ip < 0${\mathbf{ip}}<0$, or npar ≥ 20${\mathbf{npar}}\ge 20$, or npar has an invalid value, or ldcovr < npar$\mathit{ldcovr}<{\mathbf{npar}}$, or ldx < num$\mathit{ldx}<{\mathbf{num}}$, or dist ≠ 'N'${\mathbf{dist}}\ne \text{'N'}$, or dist ≠ 'T'${\mathbf{dist}}\ne \text{'T'}$, or maxit ≤ 0${\mathbf{maxit}}\le 0$, or num < max (ip,iq)${\mathbf{num}}<\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{ip}},{\mathbf{iq}}\right)$, or ${\mathbf{num}}<{\mathbf{nreg}}+{\mathbf{mn}}$.
ifail = 2${\mathbf{ifail}}=2$
 On entry, lwork < (nreg + 3) × num + 3$\mathit{lwork}<\left({\mathbf{nreg}}+3\right)×{\mathbf{num}}+3$.
ifail = 3${\mathbf{ifail}}=3$
The matrix X$X$ is not full rank.
ifail = 4${\mathbf{ifail}}=4$
The information matrix is not positive definite.
ifail = 5${\mathbf{ifail}}=5$
The maximum number of iterations has been reached.
ifail = 6${\mathbf{ifail}}=6$
The log-likelihood cannot be optimized any further.
ifail = 7${\mathbf{ifail}}=7$
No feasible model parameters could be found.

## Accuracy

Not applicable.

None.

## Example

```function nag_tsa_uni_garch_exp_estim_example
num = 100;
mn = 1;
nreg = 2;

% The series
yt = [7.53; 6.64; 7.39; 7.15; 6.42; 6.32; 6.98; 7.09; 6.63; 6.93;
7.01; 5.30; 7.86; 6.73; 7.39; 5.61; 7.02; 6.04; 7.46; 4.33;
6.02; 6.37; 3.93; 7.24; 8.58; 5.70; 9.13; 7.99; 7.79; 6.13;
8.78; 6.52; 6.79; 7.77; 7.31; 7.58; 8.78; 7.39; 8.00; 7.07;
7.65; 9.15; 8.32; 7.32; 7.58; 9.78; 8.17; 9.26; 7.79; 7.03;
7.45; 7.09; 8.06; 7.06; 9.91; 7.01; 8.32; 6.41; 8.59; 8.55;
7.77; 8.04; 9.54; 8.28; 7.97; 8.42; 8.30; 7.98; 7.60; 8.77;
7.54; 7.40; 9.26; 7.30; 9.33; 9.54; 8.08; 6.93; 4.27; 2.65;
5.03; 0.91;12.63;10.87; 9.26; 8.30; 6.85; 7.48; 9.67; 9.54;
7.33; 8.84; 7.75; 8.12; 7.29; 8.58; 7.80; 3.07; 9.33;16.91];

% The exogenous variables
x = [2.40, 0.12; 2.40, 0.12; 2.40, 0.13; 2.40, 0.14;
2.40, 0.14; 2.40, 0.15; 2.40, 0.16; 2.40, 0.16;
2.40, 0.17; 2.41, 0.18; 2.41, 0.19; 2.41, 0.19;
2.41, 0.20; 2.41, 0.21; 2.41, 0.21; 2.41, 0.22;
2.41, 0.23; 2.41, 0.23; 2.41, 0.24; 2.42, 0.25;
2.42, 0.25; 2.42, 0.26; 2.42, 0.26; 2.42, 0.27;
2.42, 0.28; 2.42, 0.28; 2.42, 0.29; 2.42, 0.30;
2.42, 0.30; 2.43, 0.31; 2.43, 0.32; 2.43, 0.32;
2.43, 0.33; 2.43, 0.33; 2.43, 0.34; 2.43, 0.35;
2.43, 0.35; 2.43, 0.36; 2.43, 0.37; 2.44, 0.37;
2.44, 0.38; 2.44, 0.38; 2.44, 0.39; 2.44, 0.39;
2.44, 0.40; 2.44, 0.41; 2.44, 0.41; 2.44, 0.42;
2.44, 0.42; 2.45, 0.43; 2.45, 0.43; 2.45, 0.44;
2.45, 0.45; 2.45, 0.45; 2.45, 0.46; 2.45, 0.46;
2.45, 0.47; 2.45, 0.47; 2.45, 0.48; 2.46, 0.48;
2.46, 0.49; 2.46, 0.49; 2.46, 0.50; 2.46, 0.50;
2.46, 0.51; 2.46, 0.51; 2.46, 0.52; 2.46, 0.52;
2.46, 0.53; 2.47, 0.53; 2.47, 0.54; 2.47, 0.54;
2.47, 0.54; 2.47, 0.55; 2.47, 0.55; 2.47, 0.56;
2.47, 0.56; 2.47, 0.57; 2.47, 0.57; 2.48, 0.57;
2.48, 0.58; 2.48, 0.58; 2.48, 0.59; 2.48, 0.59;
2.48, 0.59; 2.48, 0.60; 2.48, 0.60; 2.48, 0.61;
2.48, 0.61; 2.49, 0.61; 2.49, 0.62; 2.49, 0.62;
2.49, 0.62; 2.49, 0.63; 2.49, 0.63; 2.49, 0.63;
2.49, 0.64; 2.49, 0.64; 2.49, 0.64; 2.50, 0.64];

% Details of the model to fit
dist = 't';
ip = 1;
iq = 2;

% Control parameters
copts = true;
maxit = int64(200);
tol = 0.00001;

% Initial parameter estimates
theta = [0.05; -0.15; -0.05; 0.05; 0.15; 0.35; 3.25; 1.5; 0; 0];

% Forecast horizon
nt = 4;

% Fit the GARCH model
[theta, se, sc, covar, hp, et, ht, lgf, ifail] = ...
nag_tsa_uni_garch_exp_estim(dist, yt, x, int64(ip), int64(iq), ...
int64(nreg), int64(mn), theta, 0, copts, maxit, tol, 'ifail', int64(-1));

% Calculate the volatility forecast
[fht, ifail] = nag_tsa_uni_garch_exp_forecast(int64(nt), int64(ip), ...
int64(iq), theta, ht, et);

% Output the results
fprintf('\n               Parameter        Standard\n');
fprintf('               estimates         errors\n');

% Output the coefficient alpha_0
fprintf('Alpha0 %16.2f%16.2f\n', theta(1), se(1));
l = 2;

% Output the coefficients alpha_i
for i = l:l+iq-1
fprintf('Alpha%d %16.2f%16.2f\n', i-1, theta(i), se(i));
end
l = l+iq;

% Output the coefficients psi_i
for i = l:l+iq-1
fprintf('  Psi%d %16.2f%16.2f\n', i-l+1, theta(i), se(i));
end
l = l+iq;

% Output the coefficients beta_j
fprintf('\n');
for i = l:l+ip-1
fprintf(' Beta%d %16.2f%16.2f\n', i-l+1, theta(i), se(i));
end
l = l+ip;

% Output the estimated degrees of freedom, df
if (dist == 't')
fprintf('\n    DF %16.2f%16.2f\n', theta(l), se(l));
l = l + 1;
end

% Output the estimated mean term, b_0
if (mn == 1)
fprintf('\n    B0 %16.2f%16.2f\n', theta(l), se(l));
l = l + 1;
end

% Output the estimated linear regression coefficients, b_i
for i = l:l+nreg-1
fprintf('    B%d %16.2f%16.2f\n', i-l+1, theta(i), se(i));
end

% Display the volatility forecast
fprintf('\nVolatility forecast = %12.4f\n', fht(nt));
```
```

Parameter        Standard
estimates         errors
Alpha0             0.17            0.20
Alpha1            -0.64            0.29
Alpha2            -0.44            0.25
Psi1            -0.07            0.23
Psi2             0.35            0.24

Beta1             0.42            0.17

DF             5.65            3.69

B0           129.22           50.19
B1           -51.94           21.18
B2            13.06            3.83

Volatility forecast =       1.3373

```
```function g13fg_example
num = 100;
mn = 1;
nreg = 2;

% The series
yt = [7.53; 6.64; 7.39; 7.15; 6.42; 6.32; 6.98; 7.09; 6.63; 6.93;
7.01; 5.30; 7.86; 6.73; 7.39; 5.61; 7.02; 6.04; 7.46; 4.33;
6.02; 6.37; 3.93; 7.24; 8.58; 5.70; 9.13; 7.99; 7.79; 6.13;
8.78; 6.52; 6.79; 7.77; 7.31; 7.58; 8.78; 7.39; 8.00; 7.07;
7.65; 9.15; 8.32; 7.32; 7.58; 9.78; 8.17; 9.26; 7.79; 7.03;
7.45; 7.09; 8.06; 7.06; 9.91; 7.01; 8.32; 6.41; 8.59; 8.55;
7.77; 8.04; 9.54; 8.28; 7.97; 8.42; 8.30; 7.98; 7.60; 8.77;
7.54; 7.40; 9.26; 7.30; 9.33; 9.54; 8.08; 6.93; 4.27; 2.65;
5.03; 0.91;12.63;10.87; 9.26; 8.30; 6.85; 7.48; 9.67; 9.54;
7.33; 8.84; 7.75; 8.12; 7.29; 8.58; 7.80; 3.07; 9.33;16.91];

% The exogenous variables
x = [2.40, 0.12; 2.40, 0.12; 2.40, 0.13; 2.40, 0.14;
2.40, 0.14; 2.40, 0.15; 2.40, 0.16; 2.40, 0.16;
2.40, 0.17; 2.41, 0.18; 2.41, 0.19; 2.41, 0.19;
2.41, 0.20; 2.41, 0.21; 2.41, 0.21; 2.41, 0.22;
2.41, 0.23; 2.41, 0.23; 2.41, 0.24; 2.42, 0.25;
2.42, 0.25; 2.42, 0.26; 2.42, 0.26; 2.42, 0.27;
2.42, 0.28; 2.42, 0.28; 2.42, 0.29; 2.42, 0.30;
2.42, 0.30; 2.43, 0.31; 2.43, 0.32; 2.43, 0.32;
2.43, 0.33; 2.43, 0.33; 2.43, 0.34; 2.43, 0.35;
2.43, 0.35; 2.43, 0.36; 2.43, 0.37; 2.44, 0.37;
2.44, 0.38; 2.44, 0.38; 2.44, 0.39; 2.44, 0.39;
2.44, 0.40; 2.44, 0.41; 2.44, 0.41; 2.44, 0.42;
2.44, 0.42; 2.45, 0.43; 2.45, 0.43; 2.45, 0.44;
2.45, 0.45; 2.45, 0.45; 2.45, 0.46; 2.45, 0.46;
2.45, 0.47; 2.45, 0.47; 2.45, 0.48; 2.46, 0.48;
2.46, 0.49; 2.46, 0.49; 2.46, 0.50; 2.46, 0.50;
2.46, 0.51; 2.46, 0.51; 2.46, 0.52; 2.46, 0.52;
2.46, 0.53; 2.47, 0.53; 2.47, 0.54; 2.47, 0.54;
2.47, 0.54; 2.47, 0.55; 2.47, 0.55; 2.47, 0.56;
2.47, 0.56; 2.47, 0.57; 2.47, 0.57; 2.48, 0.57;
2.48, 0.58; 2.48, 0.58; 2.48, 0.59; 2.48, 0.59;
2.48, 0.59; 2.48, 0.60; 2.48, 0.60; 2.48, 0.61;
2.48, 0.61; 2.49, 0.61; 2.49, 0.62; 2.49, 0.62;
2.49, 0.62; 2.49, 0.63; 2.49, 0.63; 2.49, 0.63;
2.49, 0.64; 2.49, 0.64; 2.49, 0.64; 2.50, 0.64];

% Details of the model to fit
dist = 't';
ip = 1;
iq = 2;

% Control parameters
copts = true;
maxit = int64(200);
tol = 0.00001;

% Initial parameter estimates
theta = [0.05; -0.15; -0.05; 0.05; 0.15; 0.35; 3.25; 1.5; 0; 0];

% Forecast horizon
nt = 4;

% Fit the GARCH model
[theta, se, sc, covar, hp, et, ht, lgf, ifail] = ...
g13fg(dist, yt, x, int64(ip), int64(iq), int64(nreg), int64(mn), theta, 0, ...
copts, maxit, tol, 'ifail', int64(-1));

% Calculate the volatility forecast
[fht, ifail] = g13fh(int64(nt), int64(ip), int64(iq), theta, ht, et);

% Output the results
fprintf('\n               Parameter        Standard\n');
fprintf('               estimates         errors\n');

% Output the coefficient alpha_0
fprintf('Alpha0 %16.2f%16.2f\n', theta(1), se(1));
l = 2;

% Output the coefficients alpha_i
for i = l:l+iq-1
fprintf('Alpha%d %16.2f%16.2f\n', i-1, theta(i), se(i));
end
l = l+iq;

% Output the coefficients psi_i
for i = l:l+iq-1
fprintf('  Psi%d %16.2f%16.2f\n', i-l+1, theta(i), se(i));
end
l = l+iq;

% Output the coefficients beta_j
fprintf('\n');
for i = l:l+ip-1
fprintf(' Beta%d %16.2f%16.2f\n', i-l+1, theta(i), se(i));
end
l = l+ip;

% Output the estimated degrees of freedom, df
if (dist == 't')
fprintf('\n    DF %16.2f%16.2f\n', theta(l), se(l));
l = l + 1;
end

% Output the estimated mean term, b_0
if (mn == 1)
fprintf('\n    B0 %16.2f%16.2f\n', theta(l), se(l));
l = l + 1;
end

% Output the estimated linear regression coefficients, b_i
for i = l:l+nreg-1
fprintf('    B%d %16.2f%16.2f\n', i-l+1, theta(i), se(i));
end

% Display the volatility forecast
fprintf('\nVolatility forecast = %12.4f\n', fht(nt));
```
```

Parameter        Standard
estimates         errors
Alpha0             0.17            0.20
Alpha1            -0.64            0.29
Alpha2            -0.44            0.25
Psi1            -0.07            0.23
Psi2             0.35            0.24

Beta1             0.42            0.17

DF             5.65            3.69

B0           129.22           50.19
B1           -51.94           21.18
B2            13.06            3.83

Volatility forecast =       1.3373

```