home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / signal / arch_fit.m next >
Text File  |  1997-02-19  |  3KB  |  109 lines
  1. ## Copyright (C) 1995, 1996, 1997  Kurt Hornik
  2. ## 
  3. ## This program is free software; you can redistribute it and/or modify
  4. ## it under the terms of the GNU General Public License as published by
  5. ## the Free Software Foundation; either version 2, or (at your option)
  6. ## any later version.
  7. ## 
  8. ## This program is distributed in the hope that it will be useful, but
  9. ## WITHOUT ANY WARRANTY; without even the implied warranty of
  10. ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11. ## General Public License for more details. 
  12. ## 
  13. ## You should have received a copy of the GNU General Public License
  14. ## along with this file.  If not, write to the Free Software Foundation,
  15. ## 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  16.  
  17. ## usage:  [a, b] = arch_fit (y, X, p [, ITER [, gamma [, a0, b0]]]) 
  18. ##
  19. ## Fits an ARCH regression model to the time series y using the scoring
  20. ## algorithm in Engle's original ARCH paper.  The model is
  21. ##    y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),
  22. ##    h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(p+1) * e(t-p)^2,
  23. ## where e(t) is N(0, h(t)), given y up to time t-1 and X up to t.
  24. ##
  25. ## If invoked as arch_fit (y, k, p) with a positive integer k, fit an
  26. ## ARCH(k,p) process, i.e., do the above with the t-th row of X given by
  27. ## [1, y(t-1), ..., y(t-k)].
  28. ##
  29. ## Optionally, one can specify the number of iterations ITER, the
  30. ## updating factor gamma, and initial values a0 and b0 for the scoring
  31. ## algorithm.
  32. ##
  33. ## The input parameters are:
  34. ##    y     ... time series (vector)
  35. ##    X     ... matrix of (ordinary) regressors or order of
  36. ##              autoregression 
  37. ##    p     ... order of the regression of the residual variance
  38.  
  39. ## Author:  KH <Kurt.Hornik@ci.tuwien.ac.at>
  40. ## Description:  Fit an ARCH regression model
  41.  
  42. function [a, b] = arch_fit (y, X, p, ITER, gamma, a0, b0)
  43.  
  44.   if ((nargin < 3) || (nargin == 6) || (nargin > 7))
  45.     usage ("arch_fit (y, X, p [, ITER [, gamma [, a0, b0]]])");
  46.   endif
  47.   
  48.   if !(is_vector (y))
  49.     error ("arch_test:  y must be a vector");
  50.   endif
  51.  
  52.   T   = length (y);
  53.   y   = reshape (y, T, 1);
  54.   [rx, cx] = size (X);
  55.   if ((rx == 1) && (cx == 1))
  56.     X = autoreg_matrix (y, X);
  57.   elseif !(rx == T)
  58.     error (["arch_test:  ", ...
  59.         "either rows (X) == length (y), or X is a scalar"]);
  60.   endif
  61.  
  62.   [T, k] = size (X);
  63.   
  64.   if (nargin == 7)
  65.     a   = a0;
  66.     b   = b0;
  67.     e   = y - X * b;
  68.   else
  69.     [b, v_b, e] = ols (y, X);
  70.     a   = [v_b, zeros (1,p)]';
  71.     if (nargin < 5)
  72.       gamma = 0.1;
  73.       if (nargin < 4)
  74.     ITER = 50;
  75.       endif
  76.     endif
  77.   endif
  78.   
  79.   esq = e.^2;
  80.   Z   = autoreg_matrix (esq, p);
  81.  
  82.   for i = 1 : ITER;
  83.     h    = Z * a;
  84.     tmp  = esq ./ h.^2 - 1 ./ h;
  85.     s    = 1 ./ h(1:T-p);
  86.     for j = 1 : p;
  87.       s = s - a(j+1) * tmp(j+1:T-p+j);
  88.     endfor
  89.     r    = 1 ./ h(1:T-p);
  90.     for j=1:p;
  91.       r = r + 2 * h(j+1:T-p+j).^2 .* esq(1:T-p);
  92.     endfor
  93.     r   = sqrt (r);
  94.     X_tilde = X(1:T-p, :) .* (r * ones (1,k));
  95.     e_tilde = e(1:T-p) .*s ./ r;
  96.     delta_b = inv (X_tilde' * X_tilde) * X_tilde' * e_tilde;
  97.     b   = b + gamma * delta_b;
  98.     e   = y - X * b;
  99.     esq = e .^ 2;
  100.     Z   = autoreg_matrix (esq, p);
  101.     h   = Z * a;
  102.     f   = esq ./ h - ones(T,1);
  103.     Z_tilde = Z ./ (h * ones (1, p+1));
  104.     delta_a = inv (Z_tilde' * Z_tilde) * Z_tilde' * f;
  105.     a   = a + gamma * delta_a;
  106.   endfor
  107.   
  108. endfunction
  109.