/ 
PDL-Stats-0.855
River stage one • 3 direct dependents • 4 total dependents
/ lib/PDL/Stats/TS.pd

NAME

PDL::Stats::TS -- basic time series functions

DESCRIPTION

The terms FUNCTIONS and METHODS are arbitrarily used to refer to methods that are threadable and methods that are NOT threadable, respectively. Plots require PDL::Graphics::Simple.

***EXPERIMENTAL!*** In particular, bad value support is spotty and may be shaky. USE WITH DISCRETION!

SYNOPSIS

use PDL::LiteF;
use PDL::Stats::TS;

my $r = $data->acf(5);

Autocorrelation function for up to lag h. If h is not specified it's set to t-1 by default.

acf does not process bad values.

usage:

pdl> $a = sequence 10

# lags 0 .. 5

pdl> p $a->acf(5)
[1 0.7 0.41212121 0.14848485 -0.078787879 -0.25757576]
EOD
);

pp_def('acvf', Pars => 'x(t); [o]v(h)', OtherPars => 'IV lag=>h;', GenericTypes => $F, Code => ' $GENERIC(x) s, s2, m, covh; s=0; s2=0; m=0; covh=0; long T, i; T = $SIZE(t); loop(t) %{ s += $x(); s2 += $x()*$x(); %} m = s/T; loop (h) %{ if (h) { covh = 0; for (i=0; i<T-h; i++) { covh += ($x(t=>i) - m) * ($x(t=>i+h) - m); } $v() = covh; } else { $v() = s2 - T * m * m; } %} ', PMCode => pp_line_numbers(__LINE__, <<'EOF'), sub PDL::acvf { my ($self, $h) = @_; $h ||= $self->dim(0) - 1; PDL::_acvf_int($self, my $v = PDL->null, $h+1); $v; } EOF Doc => <<'EOD', =for ref

Autocovariance function for up to lag h. If h is not specified it's set to t-1 by default.

acvf does not process bad values.

usage:

pdl> $a = sequence 10

# lags 0 .. 5

pdl> p $a->acvf(5)
[82.5 57.75 34 12.25 -6.5 -21.25]

# autocorrelation

pdl> p $a->acvf(5) / $a->acvf(0)
[1 0.7 0.41212121 0.14848485 -0.078787879 -0.25757576]
EOD
);

pp_def('dseason', Pars => 'x(t); indx d(); [o]xd(t)', GenericTypes => $F, HandleBad => 1, Code => ' PDL_Indx i, max = PDL_IF_BAD(,$SIZE(t))-1, min = PDL_IF_BAD(-1,0); PDL_Indx q = ($d() % 2)? ($d() - 1) / 2 : $d() / 2; /*find good min and max ind*/ loop (t) %{ PDL_IF_BAD(if ($ISBAD($x())) continue;,) if (min < 0) min = t; max = t; %} if ($d() % 2) { loop(t) %{ PDL_IF_BAD(if (t < min || t > max) { $SETBAD(xd()); continue; },) $GENERIC(x) sum = 0; PDL_IF_BAD(PDL_Indx dd = 0;,) for (i=-q; i<=q; i++) { PDL_Indx ti = (t+i < min)? min : (t+i > max)? max : t+i ; PDL_IF_BAD(if ($ISBAD($x(t=>ti))) continue; dd++;,) sum += $x(t=>ti); } PDL_IF_BAD(if (!dd) { $SETBAD(xd()); continue; },) $xd() = sum / PDL_IF_BAD(dd,$d()); %} } else { loop(t) %{ PDL_IF_BAD(if (t < min || t > max) { $SETBAD(xd()); continue; },) $GENERIC(x) sum = 0; PDL_IF_BAD(PDL_Indx dd = 0;,) for (i=-q; i<=q; i++) { PDL_Indx ti = (t+i < min)? min : (t+i > max)? max : t+i ; PDL_IF_BAD(if ($ISBAD($x(t=>ti))) continue; dd++;,) sum += (i == q || i == -q)? .5 * $x(t=>ti) : $x(t=>ti); } PDL_IF_BAD(if (!dd) { $SETBAD(xd()); continue; } dd--; if ( ($ISBAD(x(t=>t-q)) && $ISGOOD(x(t=>t+q)) ) || ($ISBAD(x(t=>t+q)) && $ISGOOD(x(t=>t-q)) ) ) dd += .5; ,) $xd() = sum / PDL_IF_BAD(dd,$d()); %} } ', Doc => 'Deseasonalize data using moving average filter the size of period d.', );

pp_def('fill_ma', Pars => 'x(t); indx q(); [o]xf(t)', GenericTypes => $F, HandleBad => 1, Code => ' $GENERIC(x) sum, xx; PDL_Indx i, n, max = $SIZE(t) - 1; loop(t) %{ PDL_IF_BAD(if ($ISBAD(x())) { n=0; sum=0; for (i=-$q(); i<=$q(); i++) { xx = (t+i < 0)? $x(t=>0) : (t+i > max)? $x(t=>max) : $x(t=>t+i) ; if ($ISGOODVAR(xx,x)) { sum += xx; n ++; } } if (n) { $xf() = sum / n; } else { $SETBAD(xf()); } continue; },) $xf() = $x(); %} ', PMCode => pp_line_numbers(__LINE__, <<'EOF'), sub PDL::fill_ma { my ($x, $q) = @_; PDL::_fill_ma_int($x, $q, my $x_filled = PDL->null); $x_filled->check_badflag; # carp "ma window too small, still has bad value" # if $x_filled->badflag; return $x_filled; } EOF Doc => <<'EOD', =for ref

Fill missing value with moving average. xf(t) = sum(x(t-q .. t-1, t+1 .. t+q)) / 2q.

fill_ma does handle bad values. Output pdl bad flag is cleared unless the specified window size q is too small and there are still bad values. EOD );

pp_def('filter_exp', Pars => 'x(t); a(); [o]xf(t)', GenericTypes => $F, Code => ' $GENERIC(x) b, m; b = 1 - $a(); loop(t) %{ if (t) { m = $a() * $x() + b * m; } else { m = $x(); } $xf() = m; %} ', Doc => 'Filter, exponential smoothing. xf(t) = a * x(t) + (1-a) * xf(t-1)', );

pp_def('filter_ma', Pars => 'x(t); indx q(); [o]xf(t)', GenericTypes => $F, Code => ' $GENERIC(x) sum; PDL_Indx i, n, max; n = 2 * $q() + 1; max = $SIZE(t) - 1; loop(t) %{ sum = 0; for (i=-$q(); i<=$q(); i++) { sum += (t+i < 0)? $x(t=>0) : (t+i > max)? $x(t=>max) : $x(t=>t+i) ; } $xf() = sum / n; %} ', Doc => 'Filter, moving average. xf(t) = sum(x(t-q .. t+q)) / (2q + 1)', );

pp_def('mae', Pars => 'a(n); b(n); [o]c()', GenericTypes => $F, HandleBad => 1, Code => ' $GENERIC(c) sum; sum = 0; PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); loop(n) %{ PDL_IF_BAD(if ($ISBAD($a()) || $ISBAD(b())) continue; N++;,) sum += fabs( $a() - $b() ); %} if (N < 1) { $SETBAD(c()); continue; } $c() = sum / N; ', Doc => 'Mean absolute error. MAE = 1/n * sum( abs(y - y_pred) )', );

pp_def('mape', Pars => 'a(n); b(n); [o]c()', GenericTypes => $F, HandleBad => 1, Code => ' $GENERIC(c) sum; sum = 0; PDL_Indx N = PDL_IF_BAD(0,$SIZE(n)); loop(n) %{ PDL_IF_BAD(if ($ISBAD($a()) || $ISBAD(b())) continue; N++;,) sum += fabs( ($a() - $b()) / $a() ); %} if (N < 1) { $SETBAD(c()); continue; } $c() = sum / N; ', Doc => 'Mean absolute percent error. MAPE = 1/n * sum(abs((y - y_pred) / y))', );

pp_def('wmape', Pars => 'a(n); b(n); [o]c()', GenericTypes => $F, HandleBad => 1, Code => ' $GENERIC(c) sum_e=0, sum=0; loop(n) %{ PDL_IF_BAD(if ($ISBAD($a()) || $ISBAD(b())) continue;,) sum_e += fabs( $a() - $b() ); sum += fabs( $a() ); %} if (!sum) { $SETBAD(c()); continue; } $c() = sum_e / sum; ', Doc => 'Weighted mean absolute percent error. avg(abs(error)) / avg(abs(data)). Much more robust compared to mape with division by zero error (cf. Schütz, W., & Kolassa, 2006).', );

pp_def('portmanteau', Pars => 'r(h); longlong t(); [o]Q()', GenericTypes => $F, Code => ' $GENERIC(r) sum; sum = 0; loop(h) %{ if (h) sum += $r()*$r() / ($t() - h); %} $Q() = $t() * ($t()+2) * sum; ', Doc => ' =for ref

Portmanteau significance test (Ljung-Box) for autocorrelations.

Usage:

  pdl> $a = sequence 10

  # acf for lags 0-5
  # lag 0 excluded from portmanteau

  pdl> p $chisq = $a->acf(5)->portmanteau( $a->nelem )
  11.1753902662994

  # get p-value from chisq distr

  pdl> use PDL::GSL::CDF
  pdl> p 1 - gsl_cdf_chisq_P( $chisq, 5 )
  0.0480112934306748
',
);

pp_def('pred_ar', Pars => 'x(p); b(p); [o]pred(t)', OtherPars => 'IV end=>t;', GenericTypes => $F, Code => ' PDL_Indx ord = $SIZE(p); $GENERIC(x) xt, xp[ord]; loop (t) %{ if (t < ord) { xp[t] = $x(p=>t); $pred() = xp[t]; } else { xt = 0; loop(p) %{ xt += xp[p] * $b(p=>ord-p-1); xp[p] = (p < ord - 1)? xp[p+1] : xt; %} $pred() = xt; } %} ', PMCode => pp_line_numbers(__LINE__, <<'EOF'), sub PDL::pred_ar { my ($x, $b, $t, $opt) = @_; my %opt = ( CONST => 1 ); if ($opt) { $opt{uc $_} = $opt->{$_} for keys %$opt; } $b = PDL->topdl($b); # allows passing simple number my $ext; if ($opt{CONST}) { my $t_ = $t - ( $x->dim(0) - $b->dim(0) + 1 ); PDL::_pred_ar_int($x->slice([-$b->dim(0)+1,-1]), $b->slice('0:-2'), $ext = PDL->null, $t_); $ext->slice([$b->dim(0)-1,-1]) += $b->slice(-1); return $x->append( $ext->slice([$b->dim(0)-1,-1]) ); } else { my $t_ = $t - ( $x->dim(0) - $b->dim(0) ); PDL::_pred_ar_int($x->slice([-$b->dim(0),-1]), $b, $ext = PDL->null, $t_); return $x->append($ext->slice([$b->dim(0),-1])); } } EOF Doc => <<'EOD', =for ref

Calculates predicted values up to period t (extend current series up to period t) for autoregressive series, with or without constant. If there is constant, it is the last element in b, as would be returned by ols or ols_t.

pred_ar does not process bad values.

CONST  => 1,

Usage:

pdl> $x = sequence 2

  # last element is constant
pdl> $b = pdl(.8, -.2, .3)

pdl> p $x->pred_ar($b, 7)
[0       1     1.1    0.74   0.492  0.3656 0.31408]

  # no constant
pdl> p $x->pred_ar($b(0:1), 7, {const=>0})
[0       1     0.8    0.44   0.192  0.0656 0.01408]
EOD
);

pp_addpm pp_line_numbers(__LINE__, <<'EOD');

season_m

Given length of season, returns seasonal mean and variance for each period (returns seasonal mean only in scalar context).

Default options (case insensitive):

START_POSITION => 0,     # series starts at this position in season
MISSING        => -999,  # internal mark for missing points in season
PLOT  => 0,              # boolean
 # see PDL::Graphics::Simple for next options
WIN   => undef,          # pass pgswin object for more plotting control
COLOR => 1,
my ($m, $ms) = $data->season_m( 24, { START_POSITION=>2 } );

plot_dseason

Plots deseasonalized data and original data points. Opens and closes default window for plotting unless a WIN object is passed in options. Returns deseasonalized data.

Default options (case insensitive):

WIN   => undef,
COLOR => 1,        # data point color

METHODS

plot_acf

Plots and returns autocorrelations for a time series.

Default options (case insensitive):

SIG  => 0.05,      # can specify .10, .05, .01, or .001
WIN  => undef,

Usage:

pdl> $a = sequence 10

pdl> p $r = $a->plot_acf(5)
[1 0.7 0.41212121 0.14848485 -0.078787879 -0.25757576]

REFERENCES

Brockwell, P.J., & Davis, R.A. (2002). Introduction to Time Series and Forecasting (2nd ed.). New York, NY: Springer.

Schütz, W., & Kolassa, S. (2006). Foresight: advantages of the MAD/Mean ratio over the MAPE. Retrieved Jan 28, 2010, from http://www.saf-ag.com/226+M5965d28cd19.html

AUTHOR

Copyright (C) 2009 Maggie J. Xiong <maggiexyz users.sourceforge.net>

All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDL distribution.