home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21eb.zip / octave / SCRIPTS.ZIP / scripts / signal / arma_rnd.m < prev    next >
Text File  |  1997-03-04  |  2KB  |  80 lines
  1. ## Copyright (C) 1995, 1996, 1997  Friedrich Leisch
  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:  arma_rnd (a, b, v, t [, n])
  18. ##
  19. ## Returns a simulation of the ARMA model
  20. ##   x(n) = a(1) * x(n-1) + ... + a(k) * x(n-k) +
  21. ##              + e(n) + b(1) * e(n-1) + ... + b(l) * e(n-l),
  22. ## where k is the length of vector a, l is the length of vector b and e
  23. ## is gaussian white noise with variance v.  The function returns a
  24. ## vector of length t.
  25. ##
  26. ## The optional parameter n gives the number of dummy x(i) used for
  27. ## initialization, i.e., a sequence of length t+n is generated and
  28. ## x(n+1:t+n) is returned.  If n is omitted, n=100 is used. 
  29.   
  30. ## Author:  FL <Friedrich.Leisch@ci.tuwien.ac.at>
  31. ## Description:  Simulate an ARMA process
  32.  
  33. function x = arma_rnd (a, b, v, t, n)
  34.  
  35.   unwind_protect
  36.     orig_listelemok = empty_list_elements_ok;
  37.     empty_list_elements_ok = "true";
  38.   
  39.     if (nargin == 4)
  40.       n = 100;
  41.     elseif (nargin == 5)
  42.       if (!is_scalar (t))
  43.           error ("arma_rnd: n must be a scalar");
  44.       endif
  45.     else
  46.       usage ("arma_rnd (a, b, v, t [, n])");
  47.     endif
  48.  
  49.     if ( (min (size (a)) > 1) || (min (size (b)) > 1) )
  50.       error ("arma_rnd:  a and b must not be matrices");
  51.     endif
  52.   
  53.     if (!is_scalar (t))
  54.       error ("arma_rnd:  t must be a scalar");
  55.     endif
  56.   
  57.     ar = length (a);
  58.     br = length (b);
  59.  
  60.     a = reshape (a, ar, 1);
  61.     b = reshape (b, br, 1);
  62.   
  63.     a = [1; -a];            # apply our notational convention
  64.     b = [1; b];
  65.     
  66.     n = min (n, ar + br);
  67.     
  68.     e = sqrt (v) * randn (t + n, 1);
  69.   
  70.     x = filter (b, a, e);
  71.     x = x(n + 1 : t + n);
  72.  
  73.   unwind_protect_cleanup
  74.     
  75.     empty_list_elements_ok = orig_listelemok;
  76.     
  77.   end_unwind_protect
  78.  
  79. endfunction
  80.