home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts.fat / signal / arch_rnd.m < prev    next >
Text File  |  1999-04-29  |  3KB  |  87 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:  y = arch_rnd (a, b, T)  
  18. ##
  19. ## Simulates an ARCH sequence y of length T with AR coefficients b and
  20. ## CH coefficients a.
  21. ## I.e., y follows the model
  22. ##     y(t) = b(1) + b(2) * y(t-1) + ... + b(lb) * y(t-lb+1) + e(t),
  23. ## where e(t), given y up to time t-1, is N(0, h(t)), with
  24. ##     h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(la) * e(t-la+1)^2.
  25.  
  26. ## Author:  KH <Kurt.Hornik@ci.tuwien.ac.at>
  27. ## Description:  Simulate an ARCH process
  28.  
  29. function y = arch_rnd (a, b, T)
  30.   
  31.   if (nargin != 3)
  32.     usage ("arch_rnd (a, b, T)");
  33.   endif
  34.   
  35.   if !( (min (size (a)) == 1) && (min (size (b)) == 1) )
  36.     error ("arch_rnd:  a and b must both be scalars or vectors");
  37.   endif
  38.   if !( is_scal (T) && (T > 0) && (rem (T, 1) == 0) )
  39.     error ("arch_rnd:  T must be a positive integer");
  40.   endif
  41.   
  42.   if !(a(1) > 0)
  43.     error ("arch_rnd:  a(1) must be positive");
  44.   endif
  45.   ## perhaps add a test for the roots of a(z) here ...
  46.   
  47.   la = length (a);
  48.   a  = reshape (a, 1, la);
  49.   if (la == 1)
  50.     a  = [a, 0];
  51.     la = la + 1;
  52.   endif
  53.   lb = length (b);
  54.   b  = reshape (b, 1, lb);
  55.   if (lb == 1)
  56.     b  = [b, 0];
  57.     lb = lb + 1;
  58.   endif
  59.   M  = max([la, lb]);
  60.   
  61.   e  = zeros (T, 1);
  62.   h  = zeros (T, 1);
  63.   y  = zeros (T, 1);
  64.   
  65.   h(1) = a(1);
  66.   e(1) = sqrt (h(1)) * randn;
  67.   y(1) = b(1) + e(1);
  68.   
  69.   for t= 2 : M;
  70.     ta   = min ([t, la]);
  71.     h(t) = a(1) + a(2:ta) * e(t-1:t-ta+1).^2;
  72.     e(t) = sqrt (h(t)) * randn;
  73.     tb   = min ([t, lb]);
  74.     y(t) = b(1) + b(2:tb) * y(t-1:t-tb+1) + e(t);
  75.   endfor
  76.   if (T > M)
  77.     for t = M+1 : T;
  78.       h(t) = a(1) + a(2:la) * e(t-1:t-la+1).^2;
  79.       e(t) = sqrt (h(t)) * randn;
  80.       y(t) = b(1) + b(2:lb) * y(t-1:t-tb+1) + e(t);
  81.     endfor
  82.   endif
  83.   
  84.   y = y(1:T);
  85.   
  86. endfunction
  87.