1######################################################################## 2## 3## Copyright (C) 1995-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {[@var{a}, @var{b}] =} arch_fit (@var{y}, @var{x}, @var{p}, @var{iter}, @var{gamma}, @var{a0}, @var{b0}) 28## Fit an ARCH regression model to the time series @var{y} using the scoring 29## algorithm in @nospell{Engle's} original ARCH paper. 30## 31## The model is 32## 33## @example 34## @group 35## y(t) = b(1) * x(t,1) + @dots{} + b(k) * x(t,k) + e(t), 36## h(t) = a(1) + a(2) * e(t-1)^2 + @dots{} + a(p+1) * e(t-p)^2 37## @end group 38## @end example 39## 40## @noindent 41## in which @math{e(t)} is @math{N(0, h(t))}, given a time-series vector 42## @var{y} up to time @math{t-1} and a matrix of (ordinary) regressors @var{x} 43## up to @math{t}. The order of the regression of the residual variance is 44## specified by @var{p}. 45## 46## If invoked as @code{arch_fit (@var{y}, @var{k}, @var{p})} with a positive 47## integer @var{k}, fit an ARCH(@var{k}, @var{p}) process, i.e., do the above 48## with the @math{t}-th row of @var{x} given by 49## 50## @example 51## [1, y(t-1), @dots{}, y(t-k)] 52## @end example 53## 54## Optionally, one can specify the number of iterations @var{iter}, the 55## updating factor @var{gamma}, and initial values @math{a0} and @math{b0} 56## for the scoring algorithm. 57## @end deftypefn 58 59function [a, b] = arch_fit (y, x, p, iter, gamma, a0, b0) 60 61 if (nargin < 3 || nargin == 6 || nargin > 7) 62 print_usage (); 63 endif 64 65 if (! (isvector (y))) 66 error ("arch_fit: Y must be a vector"); 67 endif 68 69 T = length (y); 70 y = reshape (y, T, 1); 71 [rx, cx] = size (x); 72 if ((rx == 1) && (cx == 1)) 73 x = autoreg_matrix (y, x); 74 elseif (! (rx == T)) 75 error ("arch_fit: either rows (X) == length (Y), or X is a scalar"); 76 endif 77 78 [T, k] = size (x); 79 80 if (nargin == 7) 81 a = a0; 82 b = b0; 83 e = y - x * b; 84 else 85 [b, v_b, e] = ols (y, x); 86 a = [v_b, (zeros (1, p))]'; 87 if (nargin < 5) 88 gamma = 0.1; 89 if (nargin < 4) 90 iter = 50; 91 endif 92 endif 93 endif 94 95 esq = e.^2; 96 Z = autoreg_matrix (esq, p); 97 98 for i = 1 : iter 99 h = Z * a; 100 tmp = esq ./ h.^2 - 1 ./ h; 101 s = 1 ./ h(1:T-p); 102 for j = 1 : p 103 s -= a(j+1) * tmp(j+1:T-p+j); 104 endfor 105 r = 1 ./ h(1:T-p); 106 for j = 1:p 107 r += 2 * h(j+1:T-p+j).^2 .* esq(1:T-p); 108 endfor 109 r = sqrt (r); 110 X_tilde = x(1:T-p, :) .* (r * ones (1,k)); 111 e_tilde = e(1:T-p) .*s ./ r; 112 delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde; 113 b += gamma * delta_b; 114 e = y - x * b; 115 esq = e .^ 2; 116 Z = autoreg_matrix (esq, p); 117 h = Z * a; 118 f = esq ./ h - ones (T,1); 119 Z_tilde = Z ./ (h * ones (1, p+1)); 120 delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f; 121 a += gamma * delta_a; 122 endfor 123 124endfunction 125