2014/06/02(月)Perl Data Language 統計編 #08 「系列相関とコレログラム」
データは↓の「第２章/本文/Salary.data」
http://www.tokyo-tosho.co.jp/download/DL02122.zip
系列相関の式を見るだけでは分かりづらいので、PDLで計算＆グラフ描画です。
実行にも時間がかかるようになってきましたが、計算結果をハッシュとかでメモするなどの高速化の余地はあります。
#!/usr/bin/env perl
use strict;
use warnings;
use feature qw/say/;
use PDL::Lite;
use PDL::IO::Misc ();
use PDL::Graphics::PLplot ();
use DDP filters => { -external => [ 'PDL' ] };
my $infile = 'Salary.data';
my $colnum_time = 0;
my $colnum_tokyo = 2;
my ($yyyymm, $salary) = PDL->rcols($infile, { COLSEP => "t", INCLUDE => qr/[0-9]/ }, $colnum_time, $colnum_tokyo);
my $T = $yyyymm->nelem;
my $X = pdl [ 0 .. $T - 2 ]; # 月差
my $Y = PDL->null; # 系列相関
for my $h (0 .. $T - 2)
{
my $hensa_sekiwa = 0;
my $s_h = 0;
my $s_plus_h = 0;
for my $i ( 0 .. ($T - 1) - $h )
{
my $y_bar_h = 0;
my $y_bar_plus_h = 0;
for my $j ( 0 .. ($T - 1) - $h )
{
$y_bar_h += $salary->at($j) / ($T - $h);
$y_bar_plus_h += $salary->at($j+$h) / ($T - $h);
}
$s_h += ($salary->at($i) - $y_bar_h) ** 2;
$s_plus_h += ($salary->at($i+$h) - $y_bar_plus_h) ** 2;
$hensa_sekiwa += ($salary->at($i) - $y_bar_h) * ($salary->at($i+$h) - $y_bar_plus_h);
}
$s_h = sqrt $s_h;
$s_plus_h = sqrt $s_plus_h;
my $r_h = $hensa_sekiwa / ($s_h * $s_plus_h);
$Y = $Y->append($r_h);
say $r_h;
}
my $pl = PDL::Graphics::PLplot->new(
DEV => 'pngcairo',
FILE => 'correlogram.png',
PAGESIZE => [1800, 700],
TITLE => '東京の給与データのコレログラム',
COLOR => 'BLUE',
XLAB => '月差',
XTICK => 6,
NXSUB => 3,
YLAB => '系列相関',
);
$pl->xyplot($X, $Y);
$pl->close;
統計検定のテキストによると、縦軸を系列相関、横軸を時間差にしたものをコレログラムというらしい。テキスト通り、６ヶ月差と１２ヶ月差で周期性が見られました。
ちなみに、「PDL::Stats」のメソッドでやると↓のようになりました。系列相関係数の定義式が本によって違うんだけど、とりあえず統計検定２級のテキストに従っておけばいいのだろうか。
#!/usr/bin/env perl
use strict;
use warnings;
use feature qw/say/;
use PDL::Lite;
use PDL::Stats;
use PDL::IO::Misc ();
use PDL::Graphics::PLplot ();
use DDP filters => { -external => [ 'PDL' ] };
my $infile = 'Salary.data';
my $colnum_time = 0;
my $colnum_tokyo = 2;
my ($yyyymm, $salary) = PDL->rcols($infile, { COLSEP => "t", INCLUDE => qr/[0-9]/ }, $colnum_time, $colnum_tokyo);
my $X = pdl [ 0 .. $yyyymm->nelem - 1 ]; # 月差
my $Y = $salary->acf;
my $pl = PDL::Graphics::PLplot->new(
DEV => 'pngcairo',
FILE => 'correlogram.png',
PAGESIZE => [1800, 700],
TITLE => '東京の給与データのコレログラム',
COLOR => 'BLUE',
XLAB => '月差',
XTICK => 6,
NXSUB => 3,
YLAB => '系列相関',
BOX => [ $X->minmax, -1.0, 1.0 ],
);
$pl->xyplot($X, $Y);
$pl->close;
