line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
1
|
|
|
|
|
|
|
|
2
|
|
|
|
|
|
|
# |
3
|
|
|
|
|
|
|
# GENERATED WITH PDL::PP! Don't modify! |
4
|
|
|
|
|
|
|
# |
5
|
|
|
|
|
|
|
package PDL::Stats::TS; |
6
|
|
|
|
|
|
|
|
7
|
|
|
|
|
|
|
@EXPORT_OK = qw( PDL::PP _acf PDL::PP _acvf PDL::PP diff PDL::PP inte PDL::PP dseason PDL::PP _fill_ma PDL::PP filter_exp PDL::PP filter_ma PDL::PP mae PDL::PP mape PDL::PP wmape PDL::PP portmanteau PDL::PP _pred_ar ); |
8
|
|
|
|
|
|
|
%EXPORT_TAGS = (Func=>[@EXPORT_OK]); |
9
|
|
|
|
|
|
|
|
10
|
2
|
|
|
2
|
|
14870
|
use PDL::Core; |
|
2
|
|
|
|
|
38983
|
|
|
2
|
|
|
|
|
11
|
|
11
|
2
|
|
|
2
|
|
412
|
use PDL::Exporter; |
|
2
|
|
|
|
|
2
|
|
|
2
|
|
|
|
|
7
|
|
12
|
2
|
|
|
2
|
|
32
|
use DynaLoader; |
|
2
|
|
|
|
|
2
|
|
|
2
|
|
|
|
|
78
|
|
13
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
|
15
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
|
17
|
|
|
|
|
|
|
@ISA = ( 'PDL::Exporter','DynaLoader' ); |
18
|
|
|
|
|
|
|
push @PDL::Core::PP, __PACKAGE__; |
19
|
|
|
|
|
|
|
bootstrap PDL::Stats::TS ; |
20
|
|
|
|
|
|
|
|
21
|
|
|
|
|
|
|
|
22
|
|
|
|
|
|
|
|
23
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
|
25
|
|
|
|
|
|
|
=encoding utf8 |
26
|
|
|
|
|
|
|
|
27
|
|
|
|
|
|
|
=head1 NAME |
28
|
|
|
|
|
|
|
|
29
|
|
|
|
|
|
|
PDL::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 PDL::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 PDL::LiteF; |
40
|
|
|
|
|
|
|
use PDL::NiceSlice; |
41
|
|
|
|
|
|
|
use PDL::Stats::TS; |
42
|
|
|
|
|
|
|
|
43
|
|
|
|
|
|
|
my $r = $data->acf(5); |
44
|
|
|
|
|
|
|
|
45
|
|
|
|
|
|
|
=cut |
46
|
|
|
|
|
|
|
|
47
|
2
|
|
|
2
|
|
6
|
use Carp; |
|
2
|
|
|
|
|
15
|
|
|
2
|
|
|
|
|
78
|
|
48
|
2
|
|
|
2
|
|
358
|
use PDL::LiteF; |
|
2
|
|
|
|
|
473
|
|
|
2
|
|
|
|
|
6
|
|
49
|
2
|
|
|
2
|
|
72038
|
use PDL::NiceSlice; |
|
2
|
|
|
|
|
3
|
|
|
2
|
|
|
|
|
14
|
|
50
|
2
|
|
|
2
|
|
19695
|
use PDL::Stats::Basic; |
|
2
|
|
|
|
|
3
|
|
|
2
|
|
|
|
|
13
|
|
51
|
2
|
|
|
2
|
|
771
|
use PDL::Stats::Kmeans; |
|
2
|
|
|
|
|
3
|
|
|
2
|
|
|
|
|
12
|
|
52
|
|
|
|
|
|
|
|
53
|
|
|
|
|
|
|
$PDL::onlinedoc->scan(__FILE__) if $PDL::onlinedoc; |
54
|
|
|
|
|
|
|
|
55
|
|
|
|
|
|
|
eval { |
56
|
|
|
|
|
|
|
require PDL::Graphics::PGPLOT::Window; |
57
|
|
|
|
|
|
|
PDL::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 = \&PDL::_acf; |
81
|
|
|
|
|
|
|
|
82
|
|
|
|
|
|
|
|
83
|
|
|
|
|
|
|
|
84
|
|
|
|
|
|
|
|
85
|
|
|
|
|
|
|
|
86
|
|
|
|
|
|
|
*_acvf = \&PDL::_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 = \&PDL::acf; |
117
|
|
|
|
|
|
|
sub PDL::acf { |
118
|
2
|
|
|
2
|
0
|
2775
|
my ($self, $h) = @_; |
119
|
2
|
|
33
|
|
|
6
|
$h ||= $self->dim(0) - 1; |
120
|
2
|
|
|
|
|
70
|
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 = \&PDL::acvf; |
154
|
|
|
|
|
|
|
sub PDL::acvf { |
155
|
1
|
|
|
1
|
0
|
211
|
my ($self, $h) = @_; |
156
|
1
|
|
33
|
|
|
4
|
$h ||= $self->dim(0) - 1; |
157
|
1
|
|
|
|
|
42
|
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 = \&PDL::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 = \&PDL::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 = \&PDL::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 = \&PDL::fill_ma; |
279
|
|
|
|
|
|
|
sub PDL::fill_ma { |
280
|
1
|
|
|
1
|
0
|
5422
|
my ($x, $q) = @_; |
281
|
|
|
|
|
|
|
|
282
|
1
|
|
|
|
|
27
|
my $x_filled = $x->_fill_ma($q); |
283
|
1
|
|
|
|
|
8
|
$x_filled->check_badflag; |
284
|
|
|
|
|
|
|
|
285
|
|
|
|
|
|
|
# carp "ma window too small, still has bad value" |
286
|
|
|
|
|
|
|
# if $x_filled->badflag; |
287
|
|
|
|
|
|
|
|
288
|
1
|
|
|
|
|
61
|
return $x_filled; |
289
|
|
|
|
|
|
|
} |
290
|
|
|
|
|
|
|
|
291
|
|
|
|
|
|
|
|
292
|
|
|
|
|
|
|
|
293
|
|
|
|
|
|
|
|
294
|
|
|
|
|
|
|
|
295
|
|
|
|
|
|
|
*_fill_ma = \&PDL::_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 = \&PDL::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 = \&PDL::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 = \&PDL::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 = \&PDL::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 = \&PDL::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 PDL::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 = \&PDL::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 PDL::pred_ar { |
571
|
2
|
|
|
2
|
0
|
2039
|
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
|
|
|
|
6
|
$b = pdl $b |
576
|
|
|
|
|
|
|
unless ref $b eq 'PDL'; # allows passing simple number |
577
|
|
|
|
|
|
|
|
578
|
2
|
|
|
|
|
3
|
my $ext; |
579
|
2
|
100
|
|
|
|
6
|
if ($opt{CONST}) { |
580
|
1
|
|
|
|
|
7
|
my $t_ = $t - ( $x->dim(0) - $b->dim(0) + 1 ); |
581
|
1
|
|
|
|
|
7
|
$ext = $x(-$b->dim(0)+1:-1, )->_pred_ar($b(0:-2), $t_); |
582
|
1
|
|
|
|
|
45
|
$ext($b->dim(0)-1:-1) += $b(-1); |
583
|
1
|
|
|
|
|
39
|
return $x->append( $ext( $b->dim(0)-1 : -1 ) ); |
584
|
|
|
|
|
|
|
} |
585
|
|
|
|
|
|
|
else { |
586
|
1
|
|
|
|
|
9
|
my $t_ = $t - ( $x->dim(0) - $b->dim(0) ); |
587
|
1
|
|
|
|
|
9
|
$ext = $x(-$b->dim(0):-1, )->_pred_ar($b, $t_); |
588
|
1
|
|
|
|
|
41
|
return $x->append($ext($b->dim(0) : -1)); |
589
|
|
|
|
|
|
|
} |
590
|
|
|
|
|
|
|
} |
591
|
|
|
|
|
|
|
|
592
|
|
|
|
|
|
|
|
593
|
|
|
|
|
|
|
|
594
|
|
|
|
|
|
|
|
595
|
|
|
|
|
|
|
|
596
|
|
|
|
|
|
|
*_pred_ar = \&PDL::_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 PDL::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 PDL::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 = \&PDL::season_m; |
627
|
|
|
|
|
|
|
sub PDL::season_m { |
628
|
1
|
|
|
1
|
0
|
2165
|
my ($self, $d, $opt) = @_; |
629
|
1
|
|
|
|
|
7
|
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 PDL::Graphics::PGPLOT for more info |
635
|
|
|
|
|
|
|
COLOR => 1, |
636
|
|
|
|
|
|
|
); |
637
|
1
|
|
33
|
|
|
8
|
$opt and $opt{uc $_} = $opt->{$_} for (keys %$opt); |
638
|
1
|
50
|
33
|
|
|
5
|
if ($opt{PLOT} and !$PGPLOT) { |
639
|
0
|
|
|
|
|
0
|
carp "No PDL::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
|
|
|
|
|
203
|
my @dims = $self->dims; |
647
|
1
|
|
|
|
|
17
|
$dims[0] = $n_season * $d; |
648
|
1
|
|
|
|
|
4
|
my $data = zeroes( @dims ) + $opt{MISSING}; |
649
|
|
|
|
|
|
|
|
650
|
1
|
|
|
|
|
76
|
$data($opt{START_POSITION} : $opt{START_POSITION} + $self->dim(0)-1, ) .= $self; |
651
|
1
|
|
|
|
|
25
|
$data->badflag(1); |
652
|
1
|
|
|
|
|
3
|
$data->inplace->setvaltobad( $opt{MISSING} ); |
653
|
|
|
|
|
|
|
|
654
|
1
|
|
|
|
|
30
|
my $s = sequence $d; |
655
|
1
|
|
|
|
|
55
|
$s = $s->dummy(1, $n_season)->flat; |
656
|
1
|
|
|
|
|
31
|
$s = $s->iv_cluster(); |
657
|
|
|
|
|
|
|
|
658
|
1
|
|
|
|
|
58
|
my ($m, $ms) = $data->centroid( $s ); |
659
|
|
|
|
|
|
|
|
660
|
1
|
50
|
|
|
|
6
|
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
|
|
|
|
6
|
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 PDL::Graphics::PGPLOT for detailed graphing options. |
699
|
|
|
|
|
|
|
|
700
|
|
|
|
|
|
|
=cut |
701
|
|
|
|
|
|
|
|
702
|
|
|
|
|
|
|
*plot_dseason = \&PDL::plot_dseason; |
703
|
|
|
|
|
|
|
sub PDL::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 PDL::Graphics::PGPLOT for plotting"; |
735
|
|
|
|
|
|
|
} |
736
|
|
|
|
|
|
|
|
737
|
0
|
|
|
|
|
|
return $dsea; |
738
|
|
|
|
|
|
|
} |
739
|
|
|
|
|
|
|
|
740
|
|
|
|
|
|
|
*filt_exp = \&PDL::filt_exp; |
741
|
|
|
|
|
|
|
sub PDL::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 = \&PDL::filt_ma; |
747
|
|
|
|
|
|
|
sub PDL::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 = \&PDL::dsea; |
753
|
|
|
|
|
|
|
sub PDL::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 = \&PDL::plot_season; |
759
|
|
|
|
|
|
|
sub PDL::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 = \&PDL::plot_acf; |
795
|
|
|
|
|
|
|
sub PDL::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 PDL::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 PDL distribution. |
850
|
|
|
|
|
|
|
|
851
|
|
|
|
|
|
|
=cut |
852
|
|
|
|
|
|
|
|
853
|
|
|
|
|
|
|
|
854
|
|
|
|
|
|
|
|
855
|
|
|
|
|
|
|
; |
856
|
|
|
|
|
|
|
|
857
|
|
|
|
|
|
|
|
858
|
|
|
|
|
|
|
|
859
|
|
|
|
|
|
|
# Exit with OK status |
860
|
|
|
|
|
|
|
|
861
|
|
|
|
|
|
|
1; |
862
|
|
|
|
|
|
|
|
863
|
|
|
|
|
|
|
|