File Coverage

blib/lib/PDLA/Stats/TS.pm
Criterion Covered Total %
statement 65 125 52.0
branch 6 38 15.7
condition 5 24 20.8
subroutine 13 19 68.4
pod 0 11 0.0
total 89 217 41.0


line stmt bran cond sub pod time code
1              
2             #
3             # GENERATED WITH PDLA::PP! Don't modify!
4             #
5             package PDLA::Stats::TS;
6              
7             @EXPORT_OK = qw( PDLA::PP _acf PDLA::PP _acvf PDLA::PP diff PDLA::PP inte PDLA::PP dseason PDLA::PP _fill_ma PDLA::PP filter_exp PDLA::PP filter_ma PDLA::PP mae PDLA::PP mape PDLA::PP wmape PDLA::PP portmanteau PDLA::PP _pred_ar );
8             %EXPORT_TAGS = (Func=>[@EXPORT_OK]);
9              
10 2     2   56724 use PDLA::Core;
  2         43039  
  2         17  
11 2     2   494 use PDLA::Exporter;
  2         4  
  2         9  
12 2     2   40 use DynaLoader;
  2         3  
  2         97  
13              
14              
15              
16            
17             @ISA = ( 'PDLA::Exporter','DynaLoader' );
18             push @PDLA::Core::PP, __PACKAGE__;
19             bootstrap PDLA::Stats::TS ;
20              
21              
22              
23              
24              
25             =encoding utf8
26              
27             =head1 NAME
28              
29             PDLA::Stats::TS -- basic time series functions
30              
31             =head1 DESCRIPTION
32              
33             The terms FUNCTIONS and METHODS are arbitrarily used to refer to methods that are threadable and methods that are NOT threadable, respectively. Plots require PDLA::Graphics::PGPLOT.
34              
35             ***EXPERIMENTAL!*** In particular, bad value support is spotty and may be shaky. USE WITH DISCRETION!
36              
37             =head1 SYNOPSIS
38              
39             use PDLA::LiteF;
40             use PDLA::NiceSlice;
41             use PDLA::Stats::TS;
42              
43             my $r = $data->acf(5);
44              
45             =cut
46              
47 2     2   10 use Carp;
  2         3  
  2         119  
48 2     2   346 use PDLA::LiteF;
  2         467  
  2         12  
49 2     2   84805 use PDLA::NiceSlice;
  2         5  
  2         19  
50 2     2   29809 use PDLA::Stats::Basic;
  2         5  
  2         26  
51 2     2   1020 use PDLA::Stats::Kmeans;
  2         6  
  2         14  
52              
53             $PDLA::onlinedoc->scan(__FILE__) if $PDLA::onlinedoc;
54              
55             eval {
56             require PDLA::Graphics::PGPLOT::Window;
57             PDLA::Graphics::PGPLOT::Window->import( 'pgwin' );
58             };
59             my $PGPLOT = 1 if !$@;
60              
61             my $DEV = ($^O =~ /win/i)? '/png' : '/xs';
62              
63              
64              
65              
66              
67              
68              
69             =head1 FUNCTIONS
70              
71              
72              
73             =cut
74              
75              
76              
77              
78              
79              
80             *_acf = \&PDLA::_acf;
81              
82              
83              
84              
85              
86             *_acvf = \&PDLA::_acvf;
87              
88              
89              
90              
91             =head2 acf
92              
93             =for sig
94              
95             Signature: (x(t); int h(); [o]r(h+1))
96              
97             =for ref
98              
99             Autocorrelation function for up to lag h. If h is not specified it's set to t-1 by default.
100              
101             acf does not process bad values.
102              
103             =for usage
104              
105             usage:
106              
107             perldl> $a = sequence 10
108              
109             # lags 0 .. 5
110              
111             perldl> p $a->acf(5)
112             [1 0.7 0.41212121 0.14848485 -0.078787879 -0.25757576]
113              
114             =cut
115              
116             *acf = \&PDLA::acf;
117             sub PDLA::acf {
118 2     2 0 3642 my ($self, $h) = @_;
119 2   33     8 $h ||= $self->dim(0) - 1;
120 2         61 return $self->_acf($h+1);
121             }
122              
123             =head2 acvf
124              
125             =for sig
126              
127             Signature: (x(t); int h(); [o]v(h+1))
128              
129             =for ref
130              
131             Autocovariance function for up to lag h. If h is not specified it's set to t-1 by default.
132              
133             acvf does not process bad values.
134              
135             =for usage
136              
137             usage:
138              
139             perldl> $a = sequence 10
140              
141             # lags 0 .. 5
142              
143             perldl> p $a->acvf(5)
144             [82.5 57.75 34 12.25 -6.5 -21.25]
145              
146             # autocorrelation
147            
148             perldl> p $a->acvf(5) / $a->acvf(0)
149             [1 0.7 0.41212121 0.14848485 -0.078787879 -0.25757576]
150              
151             =cut
152              
153             *acvf = \&PDLA::acvf;
154             sub PDLA::acvf {
155 1     1 0 312 my ($self, $h) = @_;
156 1   33     3 $h ||= $self->dim(0) - 1;
157 1         52 return $self->_acvf($h+1);
158             }
159              
160              
161              
162              
163              
164             =head2 diff
165              
166             =for sig
167              
168             Signature: (x(t); [o]dx(t))
169              
170              
171             =for ref
172              
173             Differencing. DX(t) = X(t) - X(t-1), DX(0) = X(0). Can be done inplace.
174              
175              
176              
177             =for bad
178              
179             diff does not process bad values.
180             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
181              
182              
183             =cut
184              
185              
186              
187              
188              
189              
190             *diff = \&PDLA::diff;
191              
192              
193              
194              
195              
196             =head2 inte
197              
198             =for sig
199              
200             Signature: (x(n); [o]ix(n))
201              
202              
203             =for ref
204              
205             Integration. Opposite of differencing. IX(t) = X(t) + X(t-1), IX(0) = X(0). Can be done inplace.
206              
207              
208              
209             =for bad
210              
211             inte does not process bad values.
212             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
213              
214              
215             =cut
216              
217              
218              
219              
220              
221              
222             *inte = \&PDLA::inte;
223              
224              
225              
226              
227              
228             =head2 dseason
229              
230             =for sig
231              
232             Signature: (x(t); indx d(); [o]xd(t))
233              
234              
235             =for ref
236              
237             Deseasonalize data using moving average filter the size of period d.
238              
239              
240            
241              
242             =for bad
243              
244             dseason processes bad values.
245             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
246              
247              
248             =cut
249              
250              
251              
252              
253              
254              
255             *dseason = \&PDLA::dseason;
256              
257              
258              
259              
260             =head2 fill_ma
261              
262             =for sig
263              
264             Signature: (x(t); int q(); [o]xf(t))
265              
266             =for ref
267              
268             Fill missing value with moving average. xf(t) = sum(x(t-q .. t-1, t+1 .. t+q)) / 2q.
269              
270             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.
271            
272             =for usage
273              
274             my $x_filled = $x->fill_ma( $q );
275              
276             =cut
277              
278             *fill_ma = \&PDLA::fill_ma;
279             sub PDLA::fill_ma {
280 1     1 0 7527 my ($x, $q) = @_;
281              
282 1         17 my $x_filled = $x->_fill_ma($q);
283 1         9 $x_filled->check_badflag;
284              
285             # carp "ma window too small, still has bad value"
286             # if $x_filled->badflag;
287              
288 1         63 return $x_filled;
289             }
290              
291              
292              
293              
294              
295             *_fill_ma = \&PDLA::_fill_ma;
296              
297              
298              
299              
300              
301             =head2 filter_exp
302              
303             =for sig
304              
305             Signature: (x(t); a(); [o]xf(t))
306              
307              
308             =for ref
309              
310             Filter, exponential smoothing. xf(t) = a * x(t) + (1-a) * xf(t-1)
311              
312             =for usage
313              
314              
315            
316              
317             =for bad
318              
319             filter_exp does not process bad values.
320             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
321              
322              
323             =cut
324              
325              
326              
327              
328              
329              
330             *filter_exp = \&PDLA::filter_exp;
331              
332              
333              
334              
335              
336             =head2 filter_ma
337              
338             =for sig
339              
340             Signature: (x(t); indx q(); [o]xf(t))
341              
342              
343             =for ref
344              
345             Filter, moving average. xf(t) = sum(x(t-q .. t+q)) / (2q + 1)
346              
347              
348            
349              
350             =for bad
351              
352             filter_ma does not process bad values.
353             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
354              
355              
356             =cut
357              
358              
359              
360              
361              
362              
363             *filter_ma = \&PDLA::filter_ma;
364              
365              
366              
367              
368              
369             =head2 mae
370              
371             =for sig
372              
373             Signature: (a(n); b(n); float+ [o]c())
374              
375              
376              
377             =for ref
378              
379             Mean absolute error. MAE = 1/n * sum( abs(y - y_pred) )
380              
381             =for usage
382              
383             Usage:
384              
385             $mae = $y->mae( $y_pred );
386              
387              
388              
389             =for bad
390              
391             mae processes bad values.
392             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
393              
394              
395             =cut
396              
397              
398              
399              
400              
401              
402             *mae = \&PDLA::mae;
403              
404              
405              
406              
407              
408             =head2 mape
409              
410             =for sig
411              
412             Signature: (a(n); b(n); float+ [o]c())
413              
414              
415              
416             =for ref
417              
418             Mean absolute percent error. MAPE = 1/n * sum(abs((y - y_pred) / y))
419              
420             =for usage
421              
422             Usage:
423              
424             $mape = $y->mape( $y_pred );
425              
426              
427              
428             =for bad
429              
430             mape processes bad values.
431             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
432              
433              
434             =cut
435              
436              
437              
438              
439              
440              
441             *mape = \&PDLA::mape;
442              
443              
444              
445              
446              
447             =head2 wmape
448              
449             =for sig
450              
451             Signature: (a(n); b(n); float+ [o]c())
452              
453              
454              
455             =for ref
456              
457             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).
458              
459             =for usage
460              
461             Usage:
462              
463             $wmape = $y->wmape( $y_pred );
464              
465              
466              
467             =for bad
468              
469             wmape processes bad values.
470             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
471              
472              
473             =cut
474              
475              
476              
477              
478              
479              
480             *wmape = \&PDLA::wmape;
481              
482              
483              
484              
485              
486             =head2 portmanteau
487              
488             =for sig
489              
490             Signature: (r(h); longlong t(); [o]Q())
491              
492              
493             =for ref
494              
495             Portmanteau significance test (Ljung-Box) for autocorrelations.
496              
497             =for usage
498              
499             Usage:
500              
501             perldl> $a = sequence 10
502              
503             # acf for lags 0-5
504             # lag 0 excluded from portmanteau
505            
506             perldl> p $chisq = $a->acf(5)->portmanteau( $a->nelem )
507             11.1753902662994
508            
509             # get p-value from chisq distr
510              
511             perldl> use PDLA::GSL::CDF
512             perldl> p 1 - gsl_cdf_chisq_P( $chisq, 5 )
513             0.0480112934306748
514              
515              
516            
517              
518             =for bad
519              
520             portmanteau does not process bad values.
521             It will set the bad-value flag of all output piddles if the flag is set for any of the input piddles.
522              
523              
524             =cut
525              
526              
527              
528              
529              
530              
531             *portmanteau = \&PDLA::portmanteau;
532              
533              
534              
535              
536             =head2 pred_ar
537              
538             =for sig
539              
540             Signature: (x(d); b(p|p+1); int t(); [o]pred(t))
541              
542             =for ref
543              
544             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.
545              
546             pred_ar does not process bad values.
547              
548             =for options
549              
550             CONST => 1,
551              
552             =for usage
553              
554             Usage:
555              
556             perldl> $x = sequence 2
557              
558             # last element is constant
559             perldl> $b = pdl(.8, -.2, .3)
560              
561             perldl> p $x->pred_ar($b, 7)
562             [0 1 1.1 0.74 0.492 0.3656 0.31408]
563            
564             # no constant
565             perldl> p $x->pred_ar($b(0:1), 7, {const=>0})
566             [0 1 0.8 0.44 0.192 0.0656 0.01408]
567              
568             =cut
569              
570             sub PDLA::pred_ar {
571 2     2 0 2806 my ($x, $b, $t, $opt) = @_;
572 2         6 my %opt = ( CONST => 1 );
573 2   33     12 $opt and $opt{uc $_} = $opt->{$_} for (keys %$opt);
574              
575 2 50       7 $b = pdl $b
576             unless ref $b eq 'PDLA'; # allows passing simple number
577              
578 2         4 my $ext;
579 2 100       5 if ($opt{CONST}) {
580 1         7 my $t_ = $t - ( $x->dim(0) - $b->dim(0) + 1 );
581 1         9 $ext = $x(-$b->dim(0)+1:-1, )->_pred_ar($b(0:-2), $t_);
582 1         71 $ext($b->dim(0)-1:-1) += $b(-1);
583 1         45 return $x->append( $ext( $b->dim(0)-1 : -1 ) );
584             }
585             else {
586 1         18 my $t_ = $t - ( $x->dim(0) - $b->dim(0) );
587 1         6 $ext = $x(-$b->dim(0):-1, )->_pred_ar($b, $t_);
588 1         28 return $x->append($ext($b->dim(0) : -1));
589             }
590             }
591              
592              
593              
594              
595              
596             *_pred_ar = \&PDLA::_pred_ar;
597              
598              
599              
600              
601             =head2 season_m
602              
603             Given length of season, returns seasonal mean and var for each period (returns seasonal mean only in scalar context).
604              
605             =for options
606              
607             Default options (case insensitive):
608              
609             START_POSITION => 0, # series starts at this position in season
610             MISSING => -999, # internal mark for missing points in season
611             PLOT => 1, # boolean
612             # see PDLA::Graphics::PGPLOT::Window for next options
613             WIN => undef, # pass pgwin object for more plotting control
614             DEV => '/xs', # open and close dev for plotting if no WIN
615             # defaults to '/png' in Windows
616             COLOR => 1,
617              
618             See PDLA::Graphics::PGPLOT for detailed graphing options.
619              
620             =for usage
621              
622             my ($m, $ms) = $data->season_m( 24, { START_POSITION=>2 } );
623              
624             =cut
625              
626             *season_m = \&PDLA::season_m;
627             sub PDLA::season_m {
628 1     1 0 2878 my ($self, $d, $opt) = @_;
629 1         8 my %opt = (
630             START_POSITION => 0, # series starts at this position in season
631             MISSING => -999, # internal mark for missing points in season
632             PLOT => 1,
633             WIN => undef, # pass pgwin object for more plotting control
634             DEV => $DEV, # see PDLA::Graphics::PGPLOT for more info
635             COLOR => 1,
636             );
637 1   33     9 $opt and $opt{uc $_} = $opt->{$_} for (keys %$opt);
638 1 50 33     7 if ($opt{PLOT} and !$PGPLOT) {
639 0         0 carp "No PDLA::Graphics::PGPLOT, no plot :(";
640 0         0 $opt{PLOT} = 0;
641             }
642              
643 1         5 my $n_season = ($self->dim(0) + $opt{START_POSITION}) / $d;
644 1         4 $n_season = pdl($n_season)->ceil->sum;
645              
646 1         222 my @dims = $self->dims;
647 1         21 $dims[0] = $n_season * $d;
648 1         5 my $data = zeroes( @dims ) + $opt{MISSING};
649              
650 1         78 $data($opt{START_POSITION} : $opt{START_POSITION} + $self->dim(0)-1, ) .= $self;
651 1         64 $data->badflag(1);
652 1         4 $data->inplace->setvaltobad( $opt{MISSING} );
653              
654 1         28 my $s = sequence $d;
655 1         79 $s = $s->dummy(1, $n_season)->flat;
656 1         42 $s = $s->iv_cluster();
657              
658 1         55 my ($m, $ms) = $data->centroid( $s );
659              
660 1 50       9 if ($opt{PLOT}) {
661 0         0 my $w = $opt{WIN};
662 0 0       0 if (!$w) {
663 0         0 $w = pgwin( Dev=>$opt{DEV} );
664 0         0 $w->env( 0, $d-1, $m->minmax,
665             {XTitle=>'period', YTitle=>'mean'} );
666             }
667 0         0 $w->points( sequence($d), $m, {COLOR=>$opt{COLOR}, PLOTLINE=>1} );
668              
669 0 0       0 if ($m->squeeze->ndims < 2) {
670             $w->errb( sequence($d), $m, sqrt( $ms / $s->sumover ),
671 0         0 {COLOR=>$opt{COLOR}} );
672             }
673             else {
674 0         0 carp "errb does not support multi dim pdl";
675             }
676             $w->close
677 0 0       0 unless $opt{WIN};
678             }
679              
680 1 50       9 return wantarray? ($m, $ms) : $m;
681             }
682              
683             =head2 plot_dseason
684              
685             =for ref
686              
687             Plots deseasonalized data and original data points. Opens and closes default window for plotting unless a pgwin object is passed in options. Returns deseasonalized data.
688              
689             =for options
690              
691             Default options (case insensitive):
692              
693             WIN => undef,
694             DEV => '/xs', # open and close dev for plotting if no WIN
695             # defaults to '/png' in Windows
696             COLOR => 1, # data point color
697              
698             See PDLA::Graphics::PGPLOT for detailed graphing options.
699              
700             =cut
701              
702             *plot_dseason = \&PDLA::plot_dseason;
703             sub PDLA::plot_dseason {
704 0     0 0   my ($self, $d, $opt) = @_;
705 0 0         !defined($d) and croak "please set season period length";
706 0           $self = $self->squeeze;
707              
708 0           my $dsea;
709 0 0         if ($PGPLOT) {
710 0           my %opt = (
711             WIN => undef,
712             DEV => $DEV,
713             COLOR => 1, # data point color
714             );
715 0   0       $opt and $opt{uc $_} = $opt->{$_} for (keys %$opt);
716              
717 0           $dsea = $self->dsea($d);
718              
719 0           my $w = $opt{WIN};
720 0 0         if (!$opt{WIN}) {
721 0           $w = pgwin( $opt{DEV} );
722 0           $w->env( 0, $self->dim(0)-1, $self->minmax,
723             {XTitle=>'T', YTitle=>'DV'} );
724             }
725              
726 0           my $missn = ushort $self->max + 1; # ushort in case precision issue
727             $w->line( sequence($self->dim(0)), $dsea->setbadtoval( $missn ),
728 0           {COLOR=>$opt{COLOR}+1, MISSING=>$missn} );
729 0           $w->points( sequence($self->dim(0)), $self, {COLOR=>$opt{COLOR}} );
730             $w->close
731 0 0         unless $opt{WIN};
732             }
733             else {
734 0           carp "Please install PDLA::Graphics::PGPLOT for plotting";
735             }
736              
737 0           return $dsea;
738             }
739              
740             *filt_exp = \&PDLA::filt_exp;
741             sub PDLA::filt_exp {
742 0     0 0   print STDERR "filt_exp() deprecated since version 0.5.0. Please use filter_exp() instead\n";
743 0           return filter_exp( @_ );
744             }
745              
746             *filt_ma = \&PDLA::filt_ma;
747             sub PDLA::filt_ma {
748 0     0 0   print STDERR "filt_ma() deprecated since version 0.5.0. Please use filter_ma() instead\n";
749 0           return filter_ma( @_ );
750             }
751              
752             *dsea = \&PDLA::dsea;
753             sub PDLA::dsea {
754 0     0 0   print STDERR "dsea() deprecated since version 0.5.0. Please use dseason() instead\n";
755 0           return dseason( @_ );
756             }
757              
758             *plot_season = \&PDLA::plot_season;
759             sub PDLA::plot_season {
760 0     0 0   print STDERR "plot_season() deprecated since version 0.5.0. Please use season_m() instead\n";
761 0           my ($self, $d, $opt) = @_;
762 0   0       $opt and $opt{uc $_} = $opt->{$_} for (keys %$opt);
763 0           $opt->{PLOT} = 1;
764 0           return $self->season_m( $d, $opt );
765             }
766              
767             =head1 METHODS
768              
769             =head2 plot_acf
770              
771             =for ref
772              
773             Plots and returns autocorrelations for a time series.
774              
775             =for options
776              
777             Default options (case insensitive):
778              
779             SIG => 0.05, # can specify .10, .05, .01, or .001
780             DEV => '/xs', # open and close dev for plotting
781             # defaults to '/png' in Windows
782              
783             =for usage
784              
785             Usage:
786              
787             perldl> $a = sequence 10
788            
789             perldl> p $r = $a->plot_acf(5)
790             [1 0.7 0.41212121 0.14848485 -0.078787879 -0.25757576]
791              
792             =cut
793              
794             *plot_acf = \&PDLA::plot_acf;
795             sub PDLA::plot_acf {
796 0 0   0 0   my $opt = pop @_
797             if ref $_[-1] eq 'HASH';
798 0           my ($self, $h) = @_;
799 0           my $r = $self->acf($h);
800            
801 0 0         if ($PGPLOT) {
802 0           my %opt = (
803             SIG => 0.05,
804             DEV => $DEV,
805             );
806 0   0       $opt and $opt{uc $_} = $opt->{$_} for (keys %$opt);
807              
808 0           my $w = pgwin( Dev=>$opt{DEV} );
809 0           $w->env(-1, $h+1, -1.05, 1.05, {XTitle=>'lag', YTitle=>'acf'});
810 0           $w->line(pdl(-1,$h+1), zeroes(2)); # x axis
811              
812             my $y_sig = ($opt{SIG} == 0.10)? 1.64485362695147
813             : ($opt{SIG} == 0.05)? 1.95996398454005
814             : ($opt{SIG} == 0.01)? 2.5758293035489
815 0 0         : ($opt{SIG} == 0.001)? 3.29052673149193
    0          
    0          
    0          
816             : 0
817             ;
818 0 0         unless ($y_sig) {
819 0           carp "SIG outside of recognized value. default to 0.05";
820 0           $y_sig = 1.95996398454005;
821             }
822              
823 0           $w->line( pdl(-1,$h+1), ones(2) * $y_sig / sqrt($self->dim(0)),
824             { LINESTYLE=>"Dashed" } );
825 0           $w->line( pdl(-1,$h+1), ones(2) * $y_sig / sqrt($self->dim(0)) * -1,
826             { LINESTYLE=>"Dashed" } );
827 0           for my $lag (0..$h) {
828 0           $w->line( ones(2)*$lag, pdl(0, $r($lag)) );
829             }
830 0           $w->close;
831             }
832             else {
833 0           carp "Please install PDLA::Graphics::PGPLOT::Window for plotting";
834             }
835              
836 0           return $r;
837             }
838              
839             =head1 REFERENCES
840              
841             Brockwell, P.J., & Davis, R.A. (2002). Introcution to Time Series and Forecasting (2nd ed.). New York, NY: Springer.
842              
843             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
844              
845             =head1 AUTHOR
846              
847             Copyright (C) 2009 Maggie J. Xiong
848              
849             All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation as described in the file COPYING in the PDLA distribution.
850              
851             =cut
852              
853              
854              
855             ;
856              
857              
858              
859             # Exit with OK status
860              
861             1;
862              
863