octave_packages/m/signal/arch_rnd.m

   1 ## Copyright (C) 1995-2012 Kurt Hornik
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {} arch_rnd (@var{a}, @var{b}, @var{t})
  21 ## Simulate an ARCH sequence of length @var{t} with AR
  22 ## coefficients @var{b} and CH coefficients @var{a}.  I.e., the result
  23 ## @math{y(t)} follows the model
  24 ## @c Set example in small font to prevent overfull line
  25 ##
  26 ## @smallexample
  27 ## y(t) = b(1) + b(2) * y(t-1) + @dots{} + b(lb) * y(t-lb+1) + e(t),
  28 ## @end smallexample
  29 ##
  30 ## @noindent
  31 ## where @math{e(t)}, given @var{y} up to time @math{t-1}, is
  32 ## @math{N(0, h(t))}, with
  33 ## @c Set example in small font to prevent overfull line
  34 ##
  35 ## @smallexample
  36 ## h(t) = a(1) + a(2) * e(t-1)^2 + @dots{} + a(la) * e(t-la+1)^2
  37 ## @end smallexample
  38 ## @end deftypefn
  39
  40 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
  41 ## Description: Simulate an ARCH process
  42
  43 function y = arch_rnd (a, b, t)
  44
  45   if (nargin != 3)
  46     print_usage ();
  47   endif
  48
  49   if (! ((min (size (a)) == 1) && (min (size (b)) == 1)))
  50     error ("arch_rnd: A and B must both be scalars or vectors");
  51   endif
  52   if (! (isscalar (t) && (t > 0) && (rem (t, 1) == 0)))
  53     error ("arch_rnd: T must be a positive integer");
  54   endif
  55
  56   if (! (a(1) > 0))
  57     error ("arch_rnd: A(1) must be positive");
  58   endif
  59   ## perhaps add a test for the roots of a(z) here ...
  60
  61   la = length (a);
  62   a  = reshape (a, 1, la);
  63   if (la == 1)
  64     a  = [a, 0];
  65     la = la + 1;
  66   endif
  67
  68   lb = length (b);
  69   b  = reshape (b, 1, lb);
  70   if (lb == 1)
  71     b  = [b, 0];
  72     lb = lb + 1;
  73   endif
  74   m  = max([la, lb]);
  75
  76   e  = zeros (t, 1);
  77   h  = zeros (t, 1);
  78   y  = zeros (t, 1);
  79
  80   h(1) = a(1);
  81   e(1) = sqrt (h(1)) * randn;
  82   y(1) = b(1) + e(1);
  83
  84   for t = 2:m
  85     ta   = min ([t, la]);
  86     h(t) = a(1) + a(2:ta) * e(t-ta+1:t-1).^2;
  87     e(t) = sqrt (h(t)) * randn;
  88     tb   = min ([t, lb]);
  89     y(t) = b(1) + b(2:tb) * y(t-tb+1:t-1) + e(t);
  90   endfor
  91
  92   if (t > m)
  93     for t = m+1:t
  94       h(t) = a(1) + a(2:la) * e(t-la+1:t-1).^2;
  95       e(t) = sqrt (h(t)) * randn;
  96       y(t) = b(1) + b(2:lb) * y(t-tb+1:t-1) + e(t);
  97     endfor
  98   endif
  99
 100   y = y(1:t);
 101
 102 endfunction