統計数理（2012）
第 60 巻 第 1 号 159–171
c 2012 統計数理研究所


特集「時空間統計解析：新たなる分野横断的展開」［原著論文］

一般化 Whittle 法による不等間隔時空間データの
分析
松田 安昌

†

（受付 2011 年 6 月 29 日；改訂 11 月 18 日；採択 12 月 22 日）

要旨

大規模な時空間データの尤度関数を評価することは非常に困難であることが多い．ガウス尤
度関数の計算において，高次元の逆行列および行列式の演算に時間がかかるためである．そこ
で，時系列の尤度関数の近似に使われる Whittle 尤度関数を不規則に位置する時空間データに一般化した，一般化 Whittle 尤度関数を提案する．この近似尤度関数は割り算のみの演算でガウス尤度関数を非常に高速に計算することができる．一般化 Whittle 尤度関数を最大化することで得られるパラメータ推定量の一致性，漸近正規性を示し，実際の推定パフォーマンスをシ
ミュレーションおよび関東地方の地価データを用いて検証する．
キーワード： スペクトル密度関数，ピリオドグラム，離散フーリエ変換，グラム・シュ
ミット直交化，対数尤度関数，定常空間過程．
1.

はじめに

時空間データとは時間と場所の変化に応じて観測されるデータのことで，主に 2 つに分類される．ひとつが不規則な時地点において観測されるデータで point referenced data とよばれる．
多数の観測地点で観測された地価や大気汚染物質の値などが代表的な例である．もう一つが適当な区分にしたがって集計して離散データ化されたもので mesh data とよばれる．県ごとに集計されたある疾病の発病率などが典型的である．
本論考では前者の point referenced data の分析法を扱う．時空間上のランダムな地点 s1 , . . . ,sn ∈ Rd において観測された Y (s1 ), . . . , Y (sn ) に対し，次の時空間回帰モデル(1.1)

Y (sj ) = b0 + b1 α1 (sj ) + · · · + bm αm (sj ) + X(sj ) + Zj , j = 1, . . . , n,に従うとする．ここで X(s), s ∈ Rd は未観測の定常過程， Zj は平均 0，分散 σ02 の独立な誤差項，α1 (s), . . . , αm (s) は地点 s における説明変数列，b0 , . . . , bm は回帰係数である．本論考では，定常過程 X(s) の相関構造を時空間データ Y (s1 ), . . . , Y (sn ) から推定することを目的とする．な
お表現を見やすくするため観測地点 sj , j = 1, . . . , n を平面上 R2 の点として記述することにする
が，以下の議論は任意の次元に一般化される．
標本サイズが比較的少ない場合では正規尤度関数による最尤法による相関構造の分析が標準的である（Mardia and Marshall, 1984）．しかし正規尤度関数は逆行列と行列式の評価を含むため，標本サイズが大きい大規模な時空間データの最尤法によるパラメータ推定は時間がかかり†

東北大学大学院 経済学研究科：〒980–8576 宮城県仙台市青葉区川内

160

統計数理

第 60 巻

第 1 号 2012

現実的には難しい．したがって，大規模な時空間データの尤度近似法は近年の大きな関心をよぶテーマとなっている． Stein et al.（2004）は大規模な時空間データの尤度関数の近似法を提案している．これは，時系列解析で使われる自己回帰モデルの考え方を空間過程にあてはめたものである．アイデアがわかりやすく近似精度も良いので，大規模な時空間データの尤度計算
法としては標準的なものと考えることができる．その他の尤度近似法として， Kaufman et al.
（2008）
は Covariance tapering とよばれる近似法を提案している．これは，時系列解析で使われる移動平均モデルの考え方を適用したもので，ある一定の距離を超える時地点間の相関を無視
し 0 とみなして尤度近似を行う．これらの方法はいずれも空間領域において尤度関数を近似する方法と考えることができる．
本論考では空間領域に代えて，周波数領域における尤度近似法を提案する．Matsuda andYajima（2009）は時系列解析における Whittle 尤度関数を不等間隔時空間データに拡張したものを提案した．ここではこれに改良を加えた新たな尤度関数を提案する．Matsuda and Yajima（2009）で提案した Whittle 尤度関数は，フーリエ展開を用いて時空間データを周波数領域に変換する．このとき，不等間隔データにおいてはフーリエ級数が直交関数系にならないため，周
波数領域に変換する際に情報損失が生じて近似精度が落ちてしまう．そこで，フーリエ級数をそのまま使うのではなく，グラム・シュミット直交化法を用いてフーリエ関数系を直交化し，離散フーリエ変換による情報損失を最小化してから，Whittle 尤度関数を計算する方法を導く．
第 2 節ではよく知られている定常過程のスペクトル表現を紹介する．これをもとに第 3 節で一般化 Whittle 尤度関数を導出する過程を説明する．第 4 節では一般化 Whittle 法による推定
量の漸近理論を証明し，第 5，6 節ではそれぞれ関東地方の地価データ，シミュレーションデー
タを使った応用を紹介し，現実の推定パフォーマンスを明らかにする．ただし，第 4 節では，
本論考で提案する推定量そのものの漸近理論を導くことには成功しておらず，積分近似を施した推定量の漸近理論を扱うことに注意を要する．
2.

定常過程のスペクトル表現

一般化 Whittle 法を導入するために，まず定常過程のスペクトル表現を紹介する．2 次元空間上の平均 0 の定常過程 X(s)，s ∈ R2 の共分散関数 γ(u) = EX(s)X(s − u)，u ∈ R2 に対し，


(2.1)
γ(u) =
f (λ)eiu λ̇ dλR2を満たすスペクトル密度関数 f (λ)，λ ∈ R2 が存在することを仮定する．このとき，適当な仮定の下で X(s) は次のスペクトル表現を持つことが知られている（Stein, 1999, p. 23）．

(2.2)
X(s) =
exp(iλ s)M (dλ),
R2

ここで M はランダム測度とよばれるもので，任意のボレル集合 ∆ ∈ R2 に対し，
EM (∆) = 0,
EM |(∆)|2 =f (λ)dλ,∆を満たし，さらに任意の排反なボレル集合 ∆1 ，∆2 に対しEM (∆1 )M (∆2 ) = 0,を満たす．特に X(s) が実数値のとき，ランダム測度は M (dλ) = M (−dλ) を満たすことが必要

一般化 Whittle 法による不等間隔時空間データの分析

で，スペクトル表現は
(2.3)



X(s) = 2R+2

cos(iλ s)M1 (dλ) + 2


R+
2

161

sin(λ s)M2 (dλ),に帰着する．ここで R2+ は {(x, y) ∈ R2 |y > 0} であり，M1 ，M2 はそれぞれランダム測度 M の実数部および虚数部である．さらに M1 ，M2 は互いに無相関な測度で E|M1 (∆)|2 = E|M2 (∆)|2 =（Brockwell and Davis, 1991, p. 163）．
1/2 ∆ f (λ)dλ を満たす
スペクトル表現は，定常過程が様々な周波数の波の和で与えられ，その振幅が互いに無相関な確率変数となることを示している．数学的には興味深い表現であるが，このままでは統計分
析の道具としては使いにくい．積分をリーマン近似してスペクトル表現を回帰モデルの形式に書き直す．
スペクトル表現（2.3）の積分区間 R2+ を有限域 D ⊂ R2+ で近似し，さらに D を ωj ，j = 1, . . . , p
を中心とする長方形の集合 {δj ，j = 1, . . . , p} に分割する．被積分関数を δj 上で定数となるような階段関数で近似すると，スペクトル表現のリーマン近似(2.4)

X̃(s) =p


φj cos(ωj s) + ψj sin(ωj s),j=1を得る．ここで係数は φj = 2M1 (δj )，ψj = 2M2 (δj ) であり，互いに無相関でEφj = Eψj = 0,
(2.5)

E|φj |2 = E|ψj |2 = 2

f (λ)dλ  2|δj |f (ωj )δjを満たす．リーマン近似表現 X̃(s) の近似誤差 X(s) − X̃(s) を ε(s) とおくと，空間定常過程を次の回帰モデルの形式(2.6)

X(s) =p

{φj cos(ωj s) + ψj sin(ωj s)} + ε(s)j=1で表現することができる．ここで ε(s) はスペクトル表現（2.3）において D の外側の積分値とリーマン近似誤差の和となる．
3.

一般化 Whittle 尤度関数における一般化 Whittle 尤度関数を導出する．ただし，観測地点 s1 , . . . , sn
空間回帰モデル（1.1）
は長方形区間 [0, A1 ] × [0, A2 ] のなかでランダムに位置しているものとする．X(s) のスペクト
ル密度関数がパラメータ θ = (θ1 , . . . , θq ) によって f (λ; θ) とモデル化されているとき，パラメー
タ θ を推定するための一般化 Whittle 尤度関数を定義する．
スペクトル表現のリーマン近似モデル（2.6）を（1.1）に代入する．ただし，リーマン近似に使う周波数をフーリエ周波数
2πj1 2πj2
ωj =
,
A1
A2
に設定し，(j1 , j2 ) は R2+ 上の格子点でフーリエ周波数が D に含まれる範囲を動くものとし，これらの格子点によるフーリエ周波数を原点からの距離が近い順番に並べなおして，ω1 , . . . , ωp
とおくことにする．X(s) のスペクトル表現を代入した回帰モデルに対し，従属変数を y =(Y (s1 ), . . . , Y (sn )) ，独立変数を αj = (αj (s1 ), . . . , αj (sn )) ，j = 1, . . . , m，cj = (cos(ωj s1 ), . . . ,

162

統計数理

第 60 巻

第 1 号 2012

cos(ωj sn )) ，sj = (sin(ωj s1 ), . . . , sin(ωj sn )) ，j = 1, . . . , p，回帰係数を b = (b0 , . . . , bm )，Φ = (φ1 ,ψ1 , . . . , φp , ψp ) ，誤差項を ε = (ε1 , . . . , εn ) , Z = (Z1 , . . . , Zn ) とおく．さらに計画行列を B =（1.1）は行列表記(1, α1 , . . . , αm ), F = (c1 , s1 , . . . , cp , sp ) とおくと，回帰モデル(3.1)

y = Bb + F Φ + ε + Zをもつ．
B の列ベクトルと F の列ベクトルを合わせて a1 , . . . , am+2p+1 とおき，左から順にグラム・シュミットの直交化法によって正規直交基底 e1 , . . . , em+2p+1 に変換する．すなわち，
a1
,
e1 =
a1 


(3.2)
ai − i−1
j=1 (ai ej )ejei =, i = 2, . . . , m + 2p + 1,
i−1 
ai − j=1 (ai ej )ej 
によって直交化する．ここで  ·  は通常のユークリッドノルムである．最初の m + 1 個の正規直交基底をのぞいて K = (em+2 , . . . , em+2p+1 ) とおき，K  を左から
（3.1）にかけると右辺第一項が消えて，(3.3)

K  y = K  F Φ + K  (ε + Z),を得る．そして
Ψ = K yを有限フーリエ変換と定義することにする．ここで ε は，スペクトル表現（2.2）を（2.4）で近似したときの近似誤差であるが，これを平均 0，分散
τ (θ) =f (λ; θ)dλDcの独立系列で近似することにすると，（3.3）より有限フーリエ変換 Ψ = K  y の正規対数尤度関
数は，
1
1
(3.4)
log L(θ, σ02 ) = − Ψ V (θ)−1 Ψ − log V (θ)2
2となる．ここでV (θ) = K  F G(θ)F  K + (σ02 + τ (θ))I2p ,さらに，G(θ) は Φ の分散で，（2.5）よりこれを
8π 2diag (f (ω1 ; θ), f (ω1 ; θ), . . . , f (ωp ; θ), f (ωp ; θ))
A1 A2

によって近似している．
しかしながら，（3.4）を計算することは p が大きい場合，難しい．K  F が上三角行列になって V (θ) が対角行列にならないためである．そこで，本論考では K  F をその対角成分をそのままに非対角成分を 0 に置き換えた行列
Q = diagonal(K  F )で近似することを提案する．そのとき，V (θ) は対角行列になって
(3.5)

Ij =

A1 A2 Ψ2j,
16π 2 Q2jj

一般化 Whittle 法による不等間隔時空間データの分析κj (θ) =A1 A2 σ02 + τ (θ)
,
16π 2
Q2jj

163

j = 1, . . . , 2p,とおくと，（3.4）は，パラメータ θ に依存しない定数項を除いて，
p 
I2j−1
I2j
1
(3.6)
+ 1
−
2 j=1 12 f (ωj ; θ) + κ2j−1 (θ)
f
(ω
;
θ)
+ κ2j (θ)
j
2


1
1
+ log
f (ωj ; θ) + κ2j−1 (θ) + log
f (ωj ; θ) + κ2j (θ)
22に帰着する．これを一般化 Whittle 尤度関数と定義する．メッシュデータの場合のアナロジー
により I2j−1 + I2j を ωj におけるピリオドグラム I(ωj ) と定義する．なお，メッシュデータに対しては，Q2jj = n/2 となるため，一般化 Whittle 尤度関数は既存の Whittle 尤度関数に一致す
る．一般化 Whittle 尤度関数はデータ数が多い場合に有効かつ高速な近似尤度計算法を与えることが期待される．
ここで Matsuda and Yajima（2009, p. 195, Deﬁnition 2）によるピリオドグラムの定義：
IM (λ) = |J(λ)|2 ,n√J(λ) = (2π)−1 A1 A2 n−1X(sk ) exp(−λ sk ),

(3.7)

k=1

および本稿が提案するピリオドグラムとの差異についてふれておく．前者のピリオドグラム
IM (ωj ) は，
（2.6）の回帰モデルにおいて独立変数
f2j−1 = (cos(ωj s1 ), . . . , cos(ωj sn )) ,f2j = (sin(ωj s1 ), . . . , sin(ωj sn )) ,それぞれに単回帰してその最小二乗推定値を二乗して足したもの，後者は f1 , . . . , f2p を直交基
底 e1 , . . . , e2p に変換してから e2j−1 , e2j に単回帰してその最小二乗推定値を二乗して足したものに相当する．したがって，メッシュデータの場合には， f1 , . . . , f2p は直交しているので両者の定義は一致する．一方，不等間隔データの場合には前者の回帰係数は一致性を失うが，後者は直交化によってその情報損失が抑えられる．つまり不等間隔データにおいて本論考で提案する
有限フーリエ変換は（3.7）に比べて，データの情報損失が少なくなるように周波数領域に変換していることになる．
最後に p の決定についてふれる． n × 2p 行列 K の列ベクトルは直交させなければならない
ため，2p + m + 1 ≤ n の範囲で選ぶ．ただし，2p + m + 1 が n に近づくと a1 , . . . , am+2p+1 の一次独立性が保証されず，グラム・シュミット直交化の手続きが不安定になってしまう．そこで
2p + m + 1 が n に比べてやや小さくなるように p をとり，直交化法が問題なく働く範囲で決め
るとよい．
4.

漸近理論

（3.6）
を最大化する推定量 θ̂ の一致性，漸近
本節では m = 0 の場合に一般化 Whittle 尤度関数
正規性を証明する． m = 0 の場合に限定する理由は（3.3）において K  F  n/2I2p の近似が正当化されてピリオドグラムが単純な形に帰着し， Matsuda and Yajima（2009）の証明法が適用可能となるからである．ただし，m = 0 を仮定してもこのままの推定量の漸近理論を構築することができない．テイパーとよばれるものを用いて有限フーリエ変換を調整する必要がある．

164

統計数理

第 60 巻

第 1 号 2012h(x1 , x2 ) を [0, 1] × [0, 1] の外側では 0 となる関数とし，H = [0,1]×[0,1] |h(x1 , x2 )|2 dx1 dx2 としてΨ̃j = H −1/2

n

i=1

Kij yi hsi,1 si,2
,
A1 A2

,

j = 1, . . . , 2pを有限フーリエ変換とする．テイパーが定数ならば元の定義に一致する．ピリオドグラム，一
般化 Whittle 尤度の定義はこの有限フーリエ変換を使って前節と同様に与えられる．一般化
Whittle 尤度関数（3.6）をテイパーで調整してもなお漸近理論を構築することには残念ながら成功していない．漸近理論を構築するためには，（3.6）を積分近似した尤度関数
 I1 (λ)I2 (λ)
1
(4.1)
+ 1
−
2 D 12 f (λ; θ) + κ1 (θ, λ)f (λ; θ) + κ2 (θ, λ)
2


1
1
f (λ; θ) + κ1 (θ, λ) + log
f (λ; θ) + κ2 (θ, λ) dλ
+ log
2
2
を最大化する推定量を対象にしなければならない．ただし，I1 (λ)，I2 (λ)，κ1 (θ, λ)，κ2 (θ, λ) は，
λ に最も近いフーリエ周波数 ωj に対し，それぞれ（3.5）における I2j−1 ，I2j ，κ2j−1 (θ)，κ2j (θ)
を形式的に拡張したものである．ここで，I2j−1 ，I2j に必要な Ψ2j−1 ，Ψ2j を計算する際，（3.2）
において分子第一項だけをそれぞれ (cos(λs1 ), . . . , cos(λsn )) ，(sin(λs1 ), . . . , sin(λsn )) に置き換えなければならない．なお，（4.1）に対してしか漸近理論を構築することに成功していないが，現実には（3.6）を最尤推定に用いることを本稿では提案する．
漸近理論を構築する前に，漸近の意味を明らかにしておく必要がある．時系列の場合には標本サイズの増大がそのまま漸近理論につながるが，空間過程においては，標本サイズとともに観測地域の広がり方も考慮に入れる必要がある．空間過程における漸近理論にはいくつかの方法があるが，本論考では mixed asymptotics とよばれる漸近法を用いる
（仮定
（1）
）
．以下に漸近理論の構築に必要な仮定
（1）
（
– 4）
をあげる．
(1)

s1 , . . . , sn を R2 平面上の観測地点とし，各点はsj = (A1 uj,1 , A2 uj,2 )

(2)

(3)

と表されるものとする．ここで uj = (uj,1 , uj,2 )，j = 1, . . . , n は独立に [0, 1] × [0, 1] 上に一様分布する確率変数である．さらに n，A1 ，A2 は k の関数で，k → ∞ のとき，n = nk → ∞，A1 = A1 (k) → ∞，A2 = A2 (k) → ∞ かつ A1 (k)/A2 (k) = O(1)，{A1 (k)A2 (k)}3 /nk → 0 を満たす．
{Xt , t ∈ R2 } はスペクトル密度関数 f (λ)，λ = (λ1 , λ2 ) ∈ R2 をもつ平均 0 の正規定常過程であり，f (λ) は有界かつ適当な定数 C に対し，



 ∂2f
 ≤ C, p, q = 1, 2,(λ)
 ∂λp ∂λqが成り立つ．
テイパー h(x1 , x2 ) は [0, 1] × [0, 1] 内にサポートをもち，1 階偏微分可能かつ 0 ≤ m1 ，
m2 ≤ 2 に対し，
∂ m1 +m2 h(x1 , x2 )
m2
1
∂xm
1 ∂x2

(4)

が存在する．
スペクトル密度関数はパラメータ θ ∈ Θ によってモデル化され，Θ の内点 θ0 を真の値として f (θ0 , λ) = f (λ) が成り立つ．f (λ; θ) は Θ × D 上で正値をとり，λ ∈ D に対し θ に一般化 Whittle 法による不等間隔時空間データの分析

165

関して二回偏微分可能である．さらに θ1 = θ2 ならば f (θ1 , λ) = f (θ2 , λ) となる λ ∈ D の
集合は正のルベーグ測度を持つ．
定理 1. （1.1）
において，m = 0 とする．仮定 1-4 の下で k → ∞ のとき，θ̂k は θ0 に確率収束
し，かつ漸近正規性
(A1 A2 )1/2 (θ̂k − θ0 ) → N 0, 2bΓ(θ0 )−1が成立する．ここで


b=

4

[0,1]2

1
Γ(θ) = 22π|h(x)| dx

 
D

 

−2
2

[0,1]2|h(x)| dx

∂ log f (λ; θ)∂θ∂ log f (λ; θ)
∂θ

,

dλ.

以下，証明に必要な補題を述べた後，本定理の証明を行う．なお，a = (a1 , . . . , an ) ，b = (b1 , . . . ,bn ) に対し，< a, b > は内積 ni=1 ai bi を表すものとする．
補題 1. {c1 , s1 , . . . , cp , sp } = {f1 , . . . , f2p } とおく．n → ∞ のとき，
(1) n2 < fi , fi > −1 = Op (n−1/2 ), i = 1, . . . , 2p，(2) n1 < fi , fj >= Op (n−1/2 ), i = j ．
証明.

fi = (cos(ωi s1 ), . . . , cos(ωi sn )) のときを例として証明する．
< fi , fi >=ncos2 (ωi sk ) =

k=1

n
n
1
cos(2ωi sk ).
+
2
2 k=1

したがって，ωi がフーリエ周波数であることにより，

 
2
−1
< fi , fi > −1 =Ecos(2ωi u)A−1
1 A2 du = 0
n[0,A1 ]×[0,A2 ]
となり，





2
−1
−1
−1
< fi , fi > −1 = nVarcos2 (2ωi u)A−1)
1 A2 du = O(nn[0,A1 ]×[0,A2 ]であることから証明される．（2）も同様．
補題 2. fˆj を fj の f1 , . . . , fj−1 上への射影とする．
< fˆj , fˆj >= Op (n−1 A1 A2 ),< fj , fj >j = 1, . . . , 2p.

証明.
j−1< fk , fj >2< fˆj , fˆj > ≤,< fj , fj > k=1 < fk , fk >< fj , fj >j = 1, . . . , 2p.

有限区間 D に含まれるフーリエ周波数の数は高々 A1 A2 のオーダーであることと補題 1 を適用することで証明される．

統計数理

166

第 60 巻

第 1 号 2012

補題 3. IM (λ) を（3.7）で示した Matsuda and Yajima（2009, Deﬁnition 2）によるピリオドグ
ラムであるとする．（3.5）式で定義するピリオドグラム I(ωj ) = I1 (ωj ) + I2 (ωj ) を任意の λ に（4.1）のように形式的に拡張した I(λ) に対し，次の結果が成立する．

|I(λ) − IM (λ)| = op (A1 A2 )−1/2 .
証明.

まず，


2 < y, fj >
Ψj= op (A1 A2 )−1 ,
−
Qjj
n

(4.2)

j = 1, . . . , 2p,を示す．f˜j を fj の f1 , . . . , fj−1 上への射影残差とする．このとき，< y, f˜j >< y, ej >Ψj< y, f˜j >
=
=
=
˜
Qjj
< ej , fj > < fj , fj > < f˜j , f˜j >であるから，
 

 

 Ψj
< y, fj >   < y, fj >< y, fj > 
2 < y, fj >   < y, f˜j >

≤
+
−
−
−
  < f˜ , f˜ >
 Qjjn< fj , fj > < f˜j , f˜j >   < f˜j , f˜j >j j < y, fj >2 < y, fj > + −< fj , fj >n
=B1 + B2 + B3

とおく．





< f˜j − fj , f˜j − fj >< y, y >< fˆj , fˆj >== Op n−1/2 (A1 A2 )1/2 ,
B1 ≤
˜
˜
˜
˜
˜
˜
< fj , fj >< fj , fj > < fj , fj >



 
 < y, fj >  < f˜j , f˜j >   < y, fj >  < fˆj , fˆj >



1−
=
= Op n−1 (A1 A2 )1/2 ,B2 = 
< fj , fj >   < f˜j , f˜j >  < fj , fj >< f˜j , f˜j >   < y, fj >   1 − 2 < fj , fj > = Op n−1/2 (A1 A2 )−1/2 .
B3 = 
< fj , fj >n仮定（1）より，B1 , B2 , B3 は高々 op (A1 A2 )−1 であることがわかる．
（4.2）より，
A1 A2  Ψ22j−1
4 < y, f2j−1 >2  A1 A2  Ψ22j
4 < y, f2j >2 |I(ωj ) − IM (ωj )| ≤
−
−
+


16π 2  Q22j−1,2j−1n2
16π 2  Q22j,2j
n2




√
 Ψ2j−1
2 < y, f2j−1 >  √2 < y, f2j−1 > A1 A2  Ψ2j−1

−
A
A
+
=
1
2


 Q2j−1,2j−1
16π 2  Q2j−1,2j−1
n
n
√




 Ψ2j
2 < y, f2j >  √2 < y, f2j > A1 A2  Ψ2j+− A1 A2  Q2j,2j +
16π 2  Q2j,2jnn
=op (A1 A2 )−1/2√< y, y >となる．これはフーリエ周波数上で評価されたものであるが，フーリエ周波数以外にも自然に拡張される．
定理 1 の証明. 補題 3 より Lemma 7（Matsuda and Yajima, 2009）は，ピリオドグラムを本稿が定義したピリオドグラムに置き換えても成立する．一致性は Theorem1（Matsuda andYajima, 2009）と同様にして，漸近正規性は Theorem 2（Matsuda and Yajima, 2009）と同様に
Lemma7 を使って証明される．

一般化 Whittle 法による不等間隔時空間データの分析

図 1.

5.

167

関東地区公示地価データの 5573 標本地点と南西部地域（左図長方形部分）1431 標本地点の拡大図．右図の曲線は電車の路線を示す．

地価データへの応用（3.6）
を地価データの分析に応用する．地価デー
本節では本稿で提案した一般化 Whittle 尤度
タは 2001 年に国土交通省によって収集された公示地価の一部で，5573 標本地点の平米単価とその緯度と経度とともに都心からの時間距離と最寄駅からの距離とともに記録されている．図 1の左図は標本地点を単位を km とする通常の直交座標に直した散布図である．全データを定常モデルを使って分析対象とするには大きすぎるので，一部分を抜き出して分析を行う．図 1 の右図で表されている 1431 地点の地価データを対象に分析をすすめる．
地価データは場所の特性に依存した構造をもつので一定の平均値をもつ定常過程の実現値とみなすことはできない．そこで空間回帰モデル（1.1）を本地価データにあてはめることにする．
α1 (s), α2 (s) をそれぞれ地点 s における都心からの距離，最寄駅からの距離とし，地価 Y (s) を
回帰モデル
(5.1)

Y (sj ) = c0 + c1 b1 (sj ) + c2 b2 (sj ) + X(sj ) + Zj ,j = 1, . . . , 1431,

で表す．ここで X(s) は定常過程で，そのスペクトル密度関数は Matérn class（Gelfand et al.,2010, Sec. 2.7 参照）
f (λ) =

φ,(α2 + ||λ||2 )ν+1

λ ∈ R2でモデル化されているとする．Matérn class の自己共分散関数はこのスペクトル密度の逆フーリエ変換によって得られ√√ν
2 ν||h||2 ν||h||σ1, h ∈ R2 ,(5.2)
γ(h) = ν−1
Kν
2
Γ(ν)
ρ
ρ
となることが知られている．自己共分散を記述するパラメータ (ν, σ1 , ρ) はそれぞれ平滑化パラ
メータ，シル，レンジとよばれ，スペクトル密度関数を記述するパラメータ (ν, α, φ) と次の関係で結ばれている．
√
να2ν σ1
2 ν
,
φ=
.
α=
ρ
π
このモデルに対し，本論考で提案した一般化 Whittle 尤度関数（3.6）を最大化することでパラメータの推定を行う．比較のため Stein et al.（2004）
により提案された推定法（以下 Stein 法）

統計数理

168

第 60 巻

第 1 号 2012

表 1.

地価データに対する一般化 Whittle 法（Matsuda），Stein et al.
（2004）の方法，Matsuda
and Yajima（2009）の方法による最尤推定値（MLE）とその標準誤差（SE）．

図 2.

ピリオドグラムとスペクトル推定値（Matsuda（破線），Stein 法（実線））．（a）は全体図，（b）は高周波部分に拡大したもの．ここで横軸は周波数である．

および一般化 Whittle 法の改良元となった Matsuda and Yajima（2009）による推定法（以下 MY法）による推定も同時に行う．Stein 法は，大規模な時空間データの尤度関数を本論考の方法の
周波数領域ではなく時空間領域で近似を行う代表的な近似法とみなされている．なお，3 つの推定量の推定の良さを比較するため，1431 地点のデータから 100 地点をランダムに選び，残りの 1331 地点を使ってパラメータの推定を行い，推定したパラメータを使って 100 地点の予測量を構成した．パラメータ推定結果と予測誤差の二乗平均を計算した結果を表 1 に紹介している．図 2（a）（，b）にピリオドグラムと推定されたスペクトル密度関数を示す．（a）
が低周波部
分，（b）が高周波部分に拡大したものである．
まず，一般化 Whittle 法による予測誤差は MY 法のそれを改善していることから，その尤度近似の精度も向上していると考えられる．次に一般化 Whittle 法と Stein 法の推定値はかなり似通っており，前者の予測誤差は Stein 法のそれと比べてやや劣っている．また推定されたス
ペクトル密度関数を比較すると，両者の違いはほとんど見られず，違いは前者による推定値が低周波部分でややスペクトルの値が Stein 法の推定値に比べて高くなっていることくらいである．一般化 Whittle 法は Stein 法とほぼ同等な推定パフォーマンスを有する尤度近似法である
ことを示している．

一般化 Whittle 法による不等間隔時空間データの分析

6.

表 2.

Case 1: ν = 0.5，σ1 = 30.5，ρ = 4.6，σ0 = 3.46 の Matérn class で発生させたシミュ
レーションデータに対する一般化 Whittle 法（Matsuda），Stein et al.（2004）の方法，
Matsuda and Yajima（2009）の方法による最尤推定値のバイアス（bias）と平均二乗誤差の平方根（RMSE）．

表 3.

Case 2: ν = 3.0，σ1 = 30.5，ρ = 4.6，σ0 = 3.46 の Matérn class で発生させたシミュ
レーションデータに対する一般化 Whittle 法（Matsuda），Stein et al.（2004）の方法，
Matsuda and Yajima（2009）の方法による最尤推定値のバイアス（bias）と平均二乗誤差の平方根（RMSE）．

169

シミュレーション

本節は地価データ分析で実証した一般化 Whittle 法の性質をシミュレーションを使ってたしかめることを目的とする．比較対象として，前節と同様に Stein 法および MY 法を考える．MY
法は漸近的には一般化 Whittle 法と同等であることを第 4 節補題 3 において示したが，有限標本において一般化 Whittle 法がその推定パフォーマンスをどの程度改善するのかを検証する．
地価データと同じ 1431 標本地点の座標をそのまま使い，地価を Matérn class の共分散関数（5.2）において Case 1: ν = 0.5，σ1 = 30.5，ρ = 4.6，σ0 = 3.46 および Case 2: ν = 3，σ1 = 30.5，
ρ = 4.6，σ0 = 3.46 として発生させた定常データに対して，一般化 Whittle 尤度（3.6），Stein 法
および MY 法を使ってパラメータ ν, σ1 , ρ, σ0 を推定する．これを 100 回繰り返し，推定精度をバイアスと平均二乗誤差によって比較したものが表 2，3 である．ここで 3 推定量の良さを比較するため，ランダムに 100 地点を選んで残りの 1331 地点を使って推定を行い 100 地点の予測を行って予測誤差の比較も紹介している．さらに各方法の推定 1 回あたりの平均実行速度を示すと，一般化 Whittle 法，Stein 法 および MY 法はそれぞれ 11sec.，65sec.，45sec. となっている．ただ，Stein 法は尤度計算に労力を使い，一般化 Whittle 法は尤度計算の前の直交化に時間がかかる．そのため一般化 Whittle 法の尤度計算そのものは MY 法とほぼ同一の計算時間となるが，ここでは前処理の直交化を含めた計算時間を記述している．したがって，最適化のための繰り返しが多い（少ない）ほど，また標本サイズが増える（減る）ほど一般化 Whittle 法の有効性は Stein 法に比べて強まる（弱まる）ことに留意する必要がある．なお，Dell Vostro 430，
Intel Core i7 CPU，860@2.80Gh の計算機においてプログラミング言語 Ox によって各手法を記述し，その最適化パッケージ MaxBFGS によって各尤度関数最大化を行っている．

170

統計数理

第 60 巻

第 1 号 2012

シミュレーション結果は，地価データ分析で得られた一般化 Whittle 法の性質をほぼ追認す
るものとなっている．すなわち一般化 Whittle 法は MY 法のものに比べて予測誤差を改善しているが，Stein 法に比べるとやや劣っている．これは，一般化 Whittle 法の尤度近似精度はStein 法とほぼ同等の水準まで改善されていることを示す．その理由として，有限標本における MY 法は第 3 節終わりで論じた単回帰近似によって無視できない情報損失を発生させる一方，三角関数系のグラム・シュミット直交化がその情報損失を抑えることが考えられる．一般化 Whittle 法は Stein 法に比べて遜色のない高速な尤度近似法であり，かつ周波数領域における尤度近似法である．空間領域における尤度近似の代表である Stein 法の代替として，一般化
Whittle 法は大規模な時空間データの分析が望まれている現況において有用な分析方法を提供することになると考えている．

参 考 文 献Brockwell, P. and Davis, R.（1991）. Time Series: Theory and Methods, Springer-Verlag, New York.
Gelfand, A., Diggle, P., Fuentes, M. and Guttop, P.（2010）. Handbook of Spatial Statistics, CRC
Press, Boca Raton, Florida.
. Covariance tapering for likelihood based estimaKaufman, C., Schervish, M. and Nychka, D.（2008）
tion in large spatial data sets, Journal of the American Statistical Association, 103, 1556–1569.
Mardia, K. V. and Marshall, R. J.（1984）. Maximum likelihood estimation of models for residual
covariance in spatial regression, Biometrika, 71, 135–146.
. Fourier analysis of irregularly spaced data on Rd, Journal of
Matsuda, Y. and Yajima, Y.（2009）
the Royal Statistical Society Series B, 71, 191–217.
. Interpolation of Spatial Data, Springer-Verlag, New York.
Stein, M.（1999）
Stein, M. L., Chi, Z. and Welty, L. J.（2004）. Approximating likelihoods for large spatial data sets,
Journal of the Royal Statistical Society Series B, 66, 275–296.

Proceedings of the Institute of Statistical Mathematics Vol. 60, No. 1, 159–171 (2012)

171

Generalized Whittle Likelihood Method for Irregularly Spaced Data
Yasumasa Matsuda
Graduate School of Economics and Management, Tohoku University

It is sometimes diﬃcult to evaluate the likelihood function of large irregularly spaced
data. This paper proposes a method that approximates it on the frequency domain,which is considered to be a generalized version of Whittle likelihood for time series tothat for irregularly spaced data. The generalized Whittle likelihood function provides acomputer eﬃcient algorithm for evaluating likelihood functions for large irregularly spaced
data. After providing the consistency and asymptotic normality of the estimators thatmaximize the generalized Whittle likelihood function, we examine empirical performancesof the Whittle likelihood method by using simulation and real data of land prices in Kanto
area.

Key words: Spectral density function, periodogram, finite Fourier transform, Gram-Schmidt orthogonalization, log likelihood function, stationary spatial process.