home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / octa21fb.zip / octave / SCRIPTS.ZIP / scripts / signal / diffpara.m < prev    next >
Text File  |  1999-03-26  |  2KB  |  84 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:  [d, D] = diffpara (X [, a [, b]])
  18. ##
  19. ## Returns the estimator d for the differencing parameter of an
  20. ## integrated time series.
  21. ##
  22. ## The frequencies from [2*pi*a/T, 2*pi*b/T] are used for the
  23. ## estimation. If b is omitted, the interval [2*pi/T, 2*pi*a/T] is used,
  24. ## if both b and a are omitted then a = 0.5 * sqrt(T) and b = 1.5 *
  25. ## sqrt(T) is used, where T is the sample size. If X is a matrix, the
  26. ## differencing parameter of every single column is estimated.
  27. ##
  28. ## D contains the estimators for all frequencies in the intervals
  29. ## described above, d is simply mean(D).
  30. ##
  31. ## Reference: Brockwell, Peter J. & Davis, Richard A. Time Series:
  32. ## Theory and Methods Springer 1987
  33.   
  34. ## Author:  FL <Friedrich.Leisch@ci.tuwien.ac.at>
  35. ## Description:  Estimate the fractional differencing parameter
  36.  
  37. function [d, D] = diffpara (X, a, b)
  38.   
  39.   if ((nargin < 1) || (nargin > 3))
  40.     usage ("[d [, D]] = diffpara (X [, a [, b]])");
  41.   else
  42.     if is_vector (X)
  43.       n = length (X);
  44.       k = 1;
  45.       X = reshape (X, n, 1);
  46.     else
  47.       [n, k] = size(X);
  48.     endif
  49.     if (nargin == 1)
  50.       a = 0.5 * sqrt (n);
  51.       b = 1.5 * sqrt (n);
  52.     elseif (nargin == 2)
  53.       b = a;
  54.       a = 1;
  55.     endif
  56.   endif
  57.     
  58.   if !(is_scalar (a) && is_scalar (b))
  59.     error ("diffpara:  a and b must be scalars");
  60.   endif
  61.     
  62.   D = zeros (b - a + 1, k);
  63.   
  64.   for l = 1:k
  65.     
  66.     w = 2 * pi * (1 : n-1) / n;
  67.     
  68.     x = 2 * log (abs( 1 - exp (-i*w)));
  69.     y = log (periodogram (X(2:n,l)));
  70.   
  71.     x = center (x);
  72.     y = center (y);
  73.   
  74.     for m = a:b
  75.       D(m-a+1) = - x(1:m) * y(1:m) / sumsq (x(1:m));
  76.     endfor
  77.     
  78.   endfor
  79.   
  80.   d = mean (D);
  81.   
  82. endfunction
  83.  
  84.