統計数理（2003）
第 51 巻 第 2 号 199–222
c 2003 統計数理研究所


特集「個票開示問題の統計理論」［原著論文］

労働力調査とローテーション・サンプリング加納 悟†（受付 2003 年 4 月 7 日；改訂 2003 年 10 月 7 日）

要

旨

本稿では，わが国の労働力調査を例に取り，ローテーション・サンプリングによって得られた調査結果を分析するための時系列モデルを提示する．それは，失業の発生メカニズムがプロ
ビットモデルで，その潜在変数が線型状態空間モデルとして記述される，非線型時系列モデルである．それゆえ近似に依らない限り通常のカルマンフィルターのアルゴリズムを用いた最尤法で推定することはできない．本稿では，モデルの推定はギブス・サンプリング（Gibbs Sampling）
を用いて行なわれる．残念ながら，わが国の労働力調査結果の集計表にはローテーション構造に関する情報が欠落しており，モデルの有効性のチェックは現実のデータではなく擬似データによらざるを得ない．最後に，いくつかの重要と思われる応用とモデルの拡張について述べ結論とする．
キーワード： 労働力調査，ローテーション・サンプリング，プロビットモデル，ギブ
ス・サンプリング，状態空間表現．
1.

はじめに

一般に官庁統計の調査結果は集計表のかたちで公表される．そこに見られる情報は地域や年齢などの様々な特性の組み合わせに対応する調査結果の集計値である．通常統計調査のサンプ
リングには調査時点ごとに層別の無作為抽出が用いられるが，なかには一度抽出された特定の個体が繰り返し調査対象となる場合もある．わが国では労働力調査や家計調査において繰り返し調査が採用されており，諸外国においてもほとんどの労働力調査において繰り返し調査が採用されている．このような標本の抽出方法はローテーション・サンプリングと呼ばれており，調査コストの観点からは明らかに望ましい調査方法である．反面，各時点の調査結果には相関が存在し，このような特殊な標本調査方法を考慮に入れた利用が必要になる．本稿では，わが国の労働力調査を例に取り，ローテーション・サンプリングによって得られた調査結果から失業率の全国平均やその変化などの情報を得る方法について検討を加える．それによって，調査結果の公表の際，調査回数の情報が追加されることにより，研究者にとってより多様な分析が可能となることを指摘する．
オーストラリアなど諸外国においては，労働力調査をはじめとし官庁統計の調査結果の統計分析が盛んである．これに対し，わが国の統計調査結果では，プライバシー保護の観点から利用には目的外申請を義務付けており，その利用が簡単ではない．その結果，統計調査方法や調査結果の利用に関する研究は諸外国と比べかなり遅れている．
本論文の目的の 1 つは，わが国の労働力調査結果を分析するための時系列モデルを提示する†

一橋大学 経済研究所：〒186–8601 東京都国立市中 2–1; kanoh@ier.hit-u.ac.jp

200

統計数理

第 51 巻

第 2 号 2003

ことにある．本論文の構成は以下のとおりである．第 2 節において，わが国の労働力調査におけるサンプリングの特徴について簡単にまとめる．続く第 3 節では，先行研究の紹介としてオーストラリアの労働力調査の時系列モデルを紹介する．第 4 節では，わが国の労働力調査において採用されている独特のサンプリング構造を考慮に入れた時系列モデルを示す．それは，失業のメカニズムをプロビットモデルで記述し，失業のミクロデータをも分析しうる非線型時
系列モデルである．例えばある地域における高齢者の失業率の推定など，とくに少数個のサン
プルデータの分析の際，非線型性の問題は重要となる．このような場合，このモデルでは近似に依らない限り通常のカルマンフィルターのアルゴリズムを用いた最尤法で推定することはできない．それゆえ，モデルの推定はギブス・サンプリング（Gibbs Sampling）を用いて行なわれる．第 5 節では擬似データの作成とそれを用いたモデルの適合を行なう．残念ながら，わが国の労働力調査結果の集計表にはローテーション構造に関する情報が欠落しており，現実のデータを用いた分析は不可能である．それゆえ，本稿では擬似データによってローテーション構造を持つサンプリングの長所・短所について一般的に論じる．第 6 節はモデルの拡張や応用例について述べ結語とする．
2.

労働力調査とローテーション構造

わが国の，労働力調査は総務省により毎月行なわれており，その目的は，「わが国の 15 歳以
上人口について，月々の就業状態・就業時間・産業・職業などの就業状況，失業・休職の状況などの実態とその変化を把握することにより，景気判断や雇用対策等の基礎資料を得ること」（労働力調査年報より）である．そこでは全国 11 地域から国勢調査の調査区を第 1 次抽出単位，
抽出された調査区内の住戸を第 2 次抽出単位とした層化 2 段抽出法により抽出し，抽出された住戸内の 15 歳以上の世帯員を対象とする．1 調査区内の世帯数は約 50 である．その結果，得られる標本の大きさは，平成 11 年段階で調査区約 2,900，調査世帯約 40,000 世帯，15 歳以上の構成員計約 100,000 人となっている．わが国の 15 歳以上の人口は約 1 億人であることから，
労働力調査における個人の抽出率は約 0.1％（0.001）とみなされる．また労働力調査の標本設計に当たっては，（1）
前月差や前年同月差などの時系列変化の精度を高める，（2）10 地域別の結果を公表する（標本抽出は北海道，東北，南関東，北関東・甲信，北陸，東海，近畿，中国，四国，九州，沖縄の 11 地域別に行なわれるが，公表時には沖縄を九州に含め 10 地域とする）
（
，3）
調査員
の負担から 1 調査区の抽出世帯は約 15 とする，などが考慮されている（総務庁（1999）を参照）．
わが国の労働力調査の標本は以下のような特徴的な抽出方法により得られている．まず標本
調査区は 4 ヶ月継続して調査が行なわれる．その調査区は翌年同月に再び調査される．また標本調査区全体のうち 1/4 は毎月他の調査区に入れ替えられる．一方，抽出された住戸は 2 ヶ月
継続して調査される．2 ヵ月後には同一標本調査区内の他の住戸に入れ替えられるが，翌年同月には再び調査される．すなわち，各住戸は計 4 回調査される．結局，以上の組み合わせにより 8 つのグループが形成され，それらが一部ずつ入れ替えられながら調査される．各時点においては調査のタイミングについて，1 年目 1 回目，1 年目 2 回目，2 年目 1 回目，2 年目 2 回目と異なる 4 タイプが，そしてそれらがある調査区の前半に行なわれるか後半に行なわれるかによって，計 8 通りの住戸が存在することになる．図 1 はこのような調査形式をまとめたものである．
このように同一の住戸を繰り返し抽出することの利点はいくつか考えられる．第 1 に，コストが小さくて済むことである．一般に，最初の調査のコストは高く，それゆえ同じ住戸を続けて調査すれば調査費用は大きく軽減される．たとえば途上国における調査では標本を得ることが難しく，一度得た標本はできうる限り続けて利用したいことが多い．第 2 に，時点間の変化

労働力調査とローテーション・サンプリング

201

図 1. わが国の労働力調査のローテーション構造．図中の数字は通算の調査回数を表している．

をより小さな誤差で推定しうることが予想される．異なる個体の時間的変化は同一の個体の時間的変化に個体としての標本誤差を加えたものとなってしまう．とはいえ，母集団の性質は絶えず変化してゆくであろうから，それを反映するには新たな標本を採る必要がある．また，最初に選ばれた住戸が転出したり，回答しなくなるなどその質も悪化してゆくであろう．
わが国の労働力調査では特定の住戸が続けて 2 ヶ月選ばれ，その後 10 ヶ月をおいて再び 2 ヶ月選択される．McLaren and Steel（1997）ではこのようなローテーションのパターンは 2-10-2（4）と表記されている．
（4）は住戸が計 4 回選択されることを表している．労働力調査で用いられているローテーションのパターンはこの他にもいろいろあり，各国ごとに異なるローテーショ
ンのパターンが採用されている．米国の CPS（Current Population Survey）
ではわが国の労働力調査と似たローテーションの構造を採用している．各住戸は 4 ヶ月間続けて選択された後，8 ヶ月間標本からはずされ，その後再び 4 ヶ月間標本に含まれる．よって 4-8-4（8）と表記されうる．この他にも 6-6-6（12）や 1-1-1（6），1-2-1（m）などがあり，英国の四半期の労働力調査は
1-2-1（5）と考えられる．いずれも 1 年離れた時点で同一住戸を調査することにより失業率の前年同月差はより少ない誤差で推定可能となるであろう．これに対しオーストラリアでは各住戸が 8 ヶ月続けて調査され，翌年の繰り返しはない．カナダでも同様に，6 ヶ月の継続調査のみが行なわれている．
ローテーション・サンプリングにより得られた標本データの利用には注意が必要である．まず，バイアスの存在である．ローテーション・サンプリングにおいては一般に第 1 回目の調査
結果が他と比べて異なる値をとることが多い．例えば，Bailar（1975）は，様々な調査におけるローテーション・バイアスの存在を紹介し，同時に米国の Current Population Survey（CPS）においてこのようなバイアスが推定に対して与える影響について分析している．バイアスが生じる原因としては調査票の誤解や調査項目に対して正しく応えることに対する抵抗感が考えられる．我が国の労働力調査の調査票においては，一時期あいまいな表現の質問が含まれていた．
このような時期にはバイアスが生じていたとしても不思議ではない．
第 2 に，観測される標本間に相関が存在することである．そのため母集団の平均やトレンドなどのパラメータの推定には複数の時点の情報を用いることが望ましい．もし標本観測値が推定すべきパラメータの 1 次式で表現されるのであれば，標本間の相関を推定することにより，より効率的な推定量を得ることが可能であろう．すなわち一般化最小二乗法
（GLS）
により最小
分散不偏推定量を得ることが原理的に可能である．ただし，現実には標本間の相関が時間的に一定であるとは考えにくい．さらにこのような方法では時点間の相関が生じるメカニズムに関してはふれず，単に安定的な相関があると仮定することになる．
以下ではデータの生成されるメカニズムを時系列モデルを用いて定式化する．とくに 4 節ではわが国の標本観測値の生成メカニズムを状態空間表現により記述することにより，ローテー
ション・サンプリングによって得られた標本から失業率を推定する方法を提案したい．

202

3.

統計数理

第 51 巻

第 2 号 2003

オーストラリアの労働力調査の時系列モデル

オーストラリアの労働力調査では以下のような標本設計がなされている．センサスにおける調査区から標本調査区を層別抽出により選び，8 つのローテーション・グループ（RG）に分ける．各調査区から一定のルールで住戸が系統的に選ばれる．各月ごとに RG 中の 1 つのグルー
プの住戸が同一グループの他の住戸に置き換えられる．残りの住戸は再び調査される．このパターンが概ね 5 年程度継続される．
このようにオーストラリアの労働力調査では各個体は 8 ヶ月継続して調査され，その後再び調査されることはない．上述の表記法によれば，8-0-0（8）と表すことができる．このようなローテーション構造においては，当然同一の個体（住戸）が調査される 8 ヶ月間の時系列データ（集計値）には高い相関が見られる．また RG 内のデータはすべて同一の調査区に属するため，やはりある程度の相関が予想される．
Bell and Carolan（1998）はこのような相関を考慮するため，オーストラリアの労働力調査結
果を時系列モデルを用いて分析し，失業率や労働力率のトレンドを推定している．以下では，時系列モデルを用いた分析の先行研究の紹介として彼らの時系列モデルの要点をまとめる．
図 2 に与えられる，8 個のグループを考える．観測値は各グループにおける集計値 ytj（j =1, 2, . . . , 8）であり，ytj は各 RG に共通のトレンド要因 βt と季節要因 St ，そしてグループ固有jの要因 et によって構成されるとする．まずトレンド要因 βt はその変化 ∆βt = (βt − βt−1 ) が比較的スムーズに変化するものとし，∆βt = ∆βt−1 + βt

あるいはβt = 2βt−1 − βt−2 + βtと表されるとする．また季節要因 St は
St = −(St−1 + St−2 + · · · + St−11 ) + Stと表現されるものとする．ここで βt ，St はそれぞれ N (0, σβ2 )，N (0, σS2 ) に従う確率変数である．また ejt は時点 t において第 j 回目の調査を受けるようなグループの集計値 ytj の標本誤差を表すものである．より具体的に，ejt は 3 つの要因に分解される．それらは，調査区要因 Btjと個体要因 Ptj とそしてランダムな要因 Utj である．調査区要因とは標本が同一調査区にあることの影響であり，個体要因とは同一個体が抽出されることの影響を表す．これらはローテー
ション・サンプリングを行なうがゆえに考慮されなければならない要因である．ejt は 3 つの要

図 2.

オーストラリアの労働力調査のローテーション構造．図中の数字は通算の調査回数を表
している．

203

労働力調査とローテーション・サンプリング

因の合計として
ejt = Btj + Ptj + Utjと表される．さらにこれら 3 つの要因のそれぞれは以下のようなメカニズムで変動するものとする．まず，Btj は


θ B
B =
θ B
j
t

B

j−1
2
2
2 1/2 jBt
t−1 + {(1 − θB )(1 − k2 )(1 − k1 )}j = 2, 3, . . . , 8

B

8
2
2
2 1/2 jBt
t−1 + {(1 − θB )(1 − k2 )(1 − k1 )}j=1と 1 次の自己回帰モデルで表現され，月次データの相関が安定的に減衰するものと想定する．
2
ここで jBt は N (0, σB) に従う確率変数である．オーストラリアの労働力調査では同一の調査
区から 5 年間抽出が行なわれるため，各グループ内において調査区の変更の影響を考えなくともよい．よって，調査区の影響を単純に相関で記述することが可能となり，Btj は同一調査区内の個体からなる集計値が持つ相関を表す要因と解釈される．また，Ptj は，θ P+ {(1 − θ )(1 − k )}P =(1 − k ) k 
j
t

P

2
P

j−1
t−1

2
2

1/2

k1 jP t

j
2 1/2
1 Pt
2

j = 2, 3, . . . , 8j=1と表される．ここで jP t は N (0, σP2 ) に従う確率変数である．同一個体が選ばれる 8 ヶ月間には相関が存在するが，隣り合う 8 ヶ月間同士には相関がない．また，Utj については，
Utj = k2 jU t

j = 1, 2, . . . , 8

2とする．ここで jU t は N (0, σU) に従う確率変数である．また，θB ，θP ，k1 ，k2 は未知のパラメータで，Bell and Carolan（1998）ではそれらの値は集計値データの自己相関からあらかじめ推定される．
以上の想定のもとで，結局観測値 ytj は以上の要因の和で表される．すなわち，

ytj = βt + St + ejt = βt + St + Btj + Ptj + Utj

である．実際の推定には βt ，St ，Btj ，Ptj を未知の状態と考えた線型状態空間表現により，カルマン・フィルターのアルゴリズムに従って行なわれる．その結果，オーストラリアの労働力調査において θB ，θP の値は 1 に近く，ejt の分散は約 60％が Ptj によって，約 10％が Btj によって占められると同論文は報告している．残り 30％はランダム要因である．
以上，Bell and Carolan（1998）では時系列モデルを用いて失業率の値そのものを複数の要因に分解することを考えている．以下では，個人の失業の決定要因をモデル化することを考える．
4.

労働力調査結果の時系列分析

本節ではわが国の労働力調査に特有のローテーション構造 2-10-2（4）に対するモデル化を試みる．Bell and Carolan（1998）のオーストラリアの労働力調査 8-0-0（8）では標本誤差の各成分が AR モデルにより記述された．その場合，相関は安定的にかつすばやく減衰することを想定している．しかし，現実のデータに見られる相関はなかなか減衰しない．むしろランダム・ウォークとして捉えたほうがよいと思われる．また，Bell and Carolan（1998）の時系列モデル
は集計値を分析の対象とするものであるが，ここでは個々の標本の特性を記述することを考える．本来，“失業” や “労働力人口” 等の特性は質的変量であり，バイナリ-データとして記述されるものである．そこで失業のメカニズムをプロビットモデルを用いて定式化する．一般に，ある特性を満たす個体に限って分析を行なうとすれば，そのデータ数は限られたものになる．

204

統計数理

第 51 巻

第 2 号 2003

例えば，農村に住む 60 歳以上の男性の失業や転職を分析しようという場合である．このように標本数が大きくない場合には，質的データ特有の分析が必要になる．このような考え方は時系列モデルを用い小地域におけるトレンド推定を行なった Pfeﬀermann et al.（1998）と同様である．
具体的に，労働力調査の標本が調査対象に該当する特性（例えば失業，就業，転職など）を保有しているか否かは，以下のメカニズムによって決まると考える．まず，図 1 のグループjj（j = 1, 2, . . . , 8）に属する個人 i が調査特性を保有しているのであれば yit= 1，そうでなけれ
j
j∗
j∗
j
ば，yit = 0 とする．またその個人 i に対し，潜在変数 yit が存在し，yit ≥ 0 のとき，yit=1とj∗なるものとする．さらに，yit
は，個人 i の個体的特徴 αjit とトレンド要因 βt ，そしてグルー
プ要因 gtj によって構成され，
j∗yit= βt + αjit + gtjと表されるものとする．ここでグループ要因 gtj とは 8 個のそれぞれのグループに特有な要因であり，各グループに属する個体が特定の調査区から続けて抽出されていることの影響や，同一個体の調査回数に応じて生ずるローテーション・バイアスを表すものである．我が国の労働力調査においては，同一調査区から 4 回続けて同じ住戸を選ぶという工夫がなされている．その場合，同一調査区に属する個体には何らかの共通性が見られるであろうから，調査区効果が4 ヶ月続けて存在することを考慮する．よってグループ要因 gtj はバイアスと調査区効果の 2
つの和と考えればよい．その結果，同じ調査区の前半の調査結果か後半の調査結果かということと，各個体（住戸）の調査回数との組み合わせの違いで 8 つのグループに分類され，それらが図 1 の各グループに対応する．
以下では，季節性要因は考慮されていない．季節性要因はオーストラリアの労働力調査と同様に定式化することは可能である．しかし現実には，まず失業率の推定値を求めた後に，曜日
効果やうるう年効果など様々な要素を考慮に入れ，X12-ARIMA 等の移動平均法によって季節調整値が計算されるようである．Bell and Carolan（1998）はローテーション構造を考慮した場
合と X-11 により従来のように季節調整を行なった場合のトレンド推定値の比較を行なっている．わが国の労働力調査の季節調整に関しては別途議論すべき重要な課題である．以下では議論の焦点をローテーション構造に絞るため，観測値として季節調整済みの系列を用いることとし，あえて季節性要因を考慮しないことにする．
j
このとき，yit= 1 となる個体が抽出される確率は
jj∗
Pr(yit= 1) = Pr(yit≥ 0) = Pr(βt + αjit + gtj ≥ 0)jと表現される．個票レベルにおいては yit
は観察可能であるが，一般に官庁統計における個票は公開されておらず，結果は集計値として公表されるのが通常である．ここでは，各グループ
 yjにおいて Njit = Ȳtj とし，Ȳtj が観察可能であるとしよう．ただし個票レベルでデータが得られる場合にも，ここでのモデルはそのまま適用可能である．また Nj は各グループ j における標本数を表し，簡単化のため時間に関しては変化しないものとする．
j∗一方，yitの構成要因である βt ，αjit ，gtj は，それぞれ以下のメカニズムに従って時系列的に変動するものと仮定する．まずトレンド要因 βt はオーストラリアの労働力調査におけると同様，

βt = 2βt−1 − βt−2 + βtとする．ここで βt は N (0, σβ2 ) に従う確率変数である．またグループ要因 gtj の構成要素のう205

労働力調査とローテーション・サンプリング

ち調査回数に応じて生ずるローテーション・バイアス brt（r = 1, . . . , 4）はランダム・ウォークbt ∼ N (0, σb2 )brt = brt−1 + bt ;に従うものとする．このとき，ローテーション・バイアスの各グループに対する影響 bjt は

 G b + G b + G b + G b
G b +G b +G b +G bb = G b + G b + G b + G

j
t

1t 1t

2t 2t

3t 3t

4t 4t

j = 1, 5

3t 1t

4t 2t

1t 3t

2t 4t

j = 2, 64t−1 b4t

j = 3, 7G3t−1 b1t + G4t−1 b2t + G1t−1 b3t + G2t−1 b4t

j = 4, 81t−1 1t2t−1 2t3t−1 3tと表される．ここで Gkt（k = 1, . . . , 4）は以下のように定義されるダミー変数ベクトル Gk の第t 成分である．

























































12

12

12

G1 = [101010101010 | 000000000000 | 101010101010 | · · · ]
12

12

12

G2 = [010101010101 | 000000000000 | 010101010101 | · · · ]
12

12

12

G3 = [000000000000 | 101010101010 | 000000000000 | · · · ]
12

12

12

G4 = [000000000000 | 010101010101 | 000000000000 | · · · ]

ここで，
らに，

4
であり，また，G10 = G20 = G30 = 0，G40 = 1 とする．さk=1 Gkt = 1（for all t）

















































































































12

12

12

G5 = [100010001000 | 000000000000 | 100010001000 | · · · ]
12

12

12

G6 = [010001000100 | 000000000000 | 010001000100 | · · · ]
12

12

12

G7 = [001000100010 | 000000000000 | 001000100010 | · · · ]
12

12

12

G8 = [000100010001 | 000000000000 | 000100010001 | · · · ]
12

12

12

G9 = [000000000000 | 100010001000 | 000000000000 | · · · ]
12

12

12

G10 = [000000000000 | 010001000100 | 000000000000 | · · · ]
12

12

12

G11 = [000000000000 | 001000100010 | 000000000000 | · · · ]
12

12

12

G12 = [000000000000 | 000100010001 | 000000000000 | · · · ]とする．Gk（k = 5, . . . , 12）は新たな調査区の開始時点を示すためのダミー変数からなるベクトルである．このとき，各グループにおける調査区効果 dst は

206

統計数理

第 51 巻

第 2 号 2003

d1t = G5t dt + G6t d1t−1 + G7t d1t−2 + G8t d1t−3 + G9t d1t−12 + G10t d1t−13 + G11t d1t−14+ G12t d1t−15
d2t = G5t−1 dt + G6t−1 d2t−1 + G7t−1 d2t−2 + G8t−1 d2t−3 + G9t−1 d2t−12 + G10t−1 d2t−13+ G11t−1 d2t−14 + G12t−1 d2t−15
d3t = G5t−2 dt + G6t−2 d3t−1 + G7t−2 d3t−2 + G8t−2 d3t−3 + G9t−2 d3t−12 + G10t−2 d3t−13+ G11t−2 d3t−14 + G12t−2 d3t−15
d4t = G5t−3 dt + G6t−3 d4t−1 + G7t−3 d4t−2 + G8t−3 d4t−3 + G9t−3 d4t−12 + G10t−3 d4t−13+ G11t−3 d4t−14 + G12t−3 d4t−15と 4 種類に分類される．ここでも 12（for all t）が成り立つ．さらに，G9,−3 = G10,−2 =k=5 Gkt = 1
G11,−1 = G12,0 = 1 で，その他は Gkt = 0（k = 5, . . . , 8 ; t = −3, −2, −1, 0）とする．ここで dt
は時点 t において開始される新たな調査区の効果を表し，dt ∼ N (0, σd2 ) とする．
これに対し，各グループに属する個体の変動 αait は
α1it = G1t it + G2t α1it−1 + G3t α1it−12 + G4t α1it−13
α2it = G3t it + G4t α2it−1 + G1t α2it−12 + G2t α2it−13
α3it = G1t−1 it + G2t−1 α3it−1 + G3t−1 α3it−12 + G4t−1 α3it−13
α4it = G3t−1 it + G4t−1 α4it−1 + G1t−1 α4it−12 + G2t−1 α4it−13とやはり 4 種類に分類される．ここで，it は時点 t におけるサンプリングによる個体抽出の変動効果を表す．各グループに対する影響は dstr と αait の組み合わせにより 8 通りに分類される．
以上をまとめ，ローテーション・サンプリングによる観測値およびその背後にある様々な要因の変動をまとめて記述する．そのため，以下のように記号を定義する．8 つのグループからの観測値，潜在変数，そして誤差項のベクトルをそれぞれȲt = [Ȳt1 , Ȳt2 , Ȳt3 , Ȳt4 , Ȳt5 , Ȳt6 , Ȳt7 , Ȳt8 ]
1∗
2∗
3∗
4∗
5∗
6∗
7∗
8∗ 
y∗it = [yit
, yit
, yit
, yit
, yit
, yit
, yit
, yit
]

vt = [v1t , v2t , v3t , v4t , v5t , v6t , v7t , v8t ] ;vt ∼ N (0, σ 2 I8 )とする．ここで，In は n 次元の単位行列である．またローテーション・バイアス，調査区効
果，個体効果をそれぞれ
bt = [b1t , b2t , b3t , b4t ]
dt = [d1t , d1t−1 , . . . , d1t−14 |d2t , d2t−1 , . . . , d2t−14 | . . . |d4t , d4t−1 , . . . , d4t−14 ]
αt = [α1t , α1t−1 , . . . , α1t−12 |α2t , α2t−1 , . . . , α2t−12 | . . . |α4t , α4t−1 , . . . , α4t−12 ]

とし，変動要因をまとめて
θt = [βt , βt−1 , bt , dt , αt ]

と表記する．また，

207

労働力調査とローテーション・サンプリング




13

Ajt =

1
0
0

0
1
0

gjt
···
0
1

0
0

0
0

0
0

0
···
0

0
0
0

1
0

0
1


0 

0 

0 

.. 

. 
0

とする．ここで Ajt（j = 1, . . . , 4）は各グループに対する個体変動の影響を表す 13×13 行列であり，その第 1 行はそれぞれ







13

g1t = [G2t , 0, . . . , 0, G3t , G4t ]











13

g2t = [G4t , 0, . . . , 0, G1t , G2t ] 





13

g3t = [G2t−1 , 0, . . . , 0, G3t−1 , G4t−1 ]









13

g4t = [G4t−1 , 0, . . . , 0, G1t−1 , G2t−1 ] となっている．同様に，




15

Djt =

1
0
0

0
1
0

hjt
···
0
1

0
0

0
0

0
0

0
···
0

0
0
0

1
0

0
1



0 

0 

0 

.. 

. 
0

とおく．ここで Djt（j = 1, . . . , 4）は調査区効果を表す 15×15 行列であり，その第 1 行はそれ
ぞれ







15

h1t = [G6t , G7t , G8t , 0, . . . , 0, G9t , G10t , G11t , G12t ] 



















15

h2t = [G6t−1 , G7t−1 , G8t−1 , 0, . . . , 0, G9t−1 , G10t−1 , G11t−1 , G12t−1 ] 
15

h3t = [G6t−2 , G7t−2 , G8t−2 , 0, . . . , 0, G9t−2 , G10t−2 , G11t−2 , G12t−2 ] 
15

h4t = [G6t−3 , G7t−3 , G8t−3 , 0, . . . , 0, G9t−3 , G10t−3 , G11t−3 , G12t−3 ] となっている．さらに表記の簡便化のため，以下のようなベクトルを定義する．













13



G1 = [G1t , 0, . . . , 0]G−
1 = [G1t−1 , 0, . . . , 0]15G5 = [G5t , 0, . . . , 0] 











13



G3 = [G3t , 0, . . . , 0] 



13





13


G−
3 = [G3t−1 , 0, . . . , 0]15



G6 = [G6t , 0, . . . , 0] 

208

統計数理





第 51 巻

15

G7 = [G7t , 0, . . . , 0]

第 2 号 2003









15

G8 = [G8t , 0, . . . , 0]     

1512

ι2 = [1, 0, · · · , 0]

ι1 = [1, 0, . . . , 0]以上の準備のもとで，ローテーション・サンプリングの時系列モデルは以下のようにまとめら
れる．
(4.1)

j
j∗
Pr(yit= 1) = Pr(yit≥ 0),Ȳtj =



j
yit
/Nj

i

(4.2)

1
1
1
1
y ∗it =
1
1
1
1

0
0
0
0
0
0
0
0

G1t
G3t
G1t−1
G3t−1
G1t
G3t
G1t−1
G3t−1

= Hθt + v t ;

G2t
G4t
G2t−1
G4t−1
G2t
G4t
G2t−1
G4t−1

G3t
G1t
G3t−1
G1t−1
G3t
G1t
G3t−1
G1t−1

G4t
G2t
G4t−1
G2t−1
G4t
G2t
G4t−1
G2t−1

ι1
ι1

ι2
0

0



ι

0
ι 
 θ + v

0

ι

ι
ι2

ι1
ι1

2

2

ι1
ι1
0

ι2

t

t

2

ι1
ι1

2

0

ι2

v t ∼ N (0, R)

1∗
2∗
3∗
4∗
5∗
6∗
7∗
8∗ 
ただし，y ∗it = [yit
, yit
, yit
, yit
, yit
, yit
, yit
, yit
] ，R = σ 2 I8 で H は 8 × 118 の行列，θt は 118∗
次元のベクトルである．ここで潜在変数 y it は状態変数の線型方程式 (4.2) として表現されるが，(4.2) よりȲt∗ = Hθt + v̄ t ,(4.2)

v̄ t ∼ N (0, R̄)が成り立つ．ここで，Ȳt∗ = [Ȳt1∗ , Ȳt2∗ , Ȳt3∗ , Ȳt4∗ , Ȳt5∗ , Ȳt6∗ , Ȳt7∗ , Ȳt8∗ ] , Ȳtj∗ =






1/N1
0

1/N2
R̄ = σ 2

..

.

0

j∗
i yit /Nj

1/N7
1/N8

である．そこで (4.1) より Ȳt∗ を発生し，(4.2) を観測方程式とみなす．また状態方程式は






 θ






2 −1
1

0

0

I4
D1t
D2t
(4.3)

θt =

D3t
D4t
A1t
A2t
0

A3t
A4t

1 0 0
0 0 0
0 1 0
0 0 G50 0 G6t−1 + 0 0 G7
0 0 G8
0 0 0
0 0 0
0 0 0
0 0 0

0
0
0
0
0
0
0
G1
G3
G−
1
G−
3




 
     
 



βt
bt

dt
it

労働力調査とローテーション・サンプリング209

βt ∼ N (0, σβ2 ), it ∼ N (0, σα2 )= Ft θt−1 + Jt ¯t ;となる．ここで，Ft は 118 × 118 の行列，Jt は 118 × 4 の行列である．
以上のシステムの推定は，観測値を得るメカニズム（4.1）が非線型であるため，通常のカル
マン・フィルターのアルゴリズムを用いた最尤法で推定することはできない．Kim and Nelson（1999）は状態空間モデルにおいて Gibbs Sampling を用いてベイズ法により状態を推定するアルゴリズムを提案している．以下の方法はこれを拡張したものである．特に以下にあるように∗（4.1）を利用して発生した Ȳt∗ を（4.2）の yitに置換えパラメータ推定を行なう．
具体的には，推定すべきパラメーター・ベクトルをΘ = [β1 , β2 , . . . , βT , b11 , . . . , b1T , . . . , b41 , . . . , b4T , d11 , . . . , d1T , . . . , d41 , . . . , d4T ,α11 , . . . , α1T , . . . , α41 , . . . , α4T , σβ2 , σb2 , σd2 , σα2 , σ 2 ] = β̃T , γ̃T , σ̃ 2 = θ̃T , σ̃ 2と設定し，以下のステップ
（a）
（
，b）
を繰り返す．
(k+1)(a) g1 (θ̃T |σ̃ 2(k) , Ỹ T ) から，θ̃Tを生成する．
(k+1)(b) g2 (σ̃ 2 |θ̃T
, Ỹ T ) から，σ̃ 2(k+1) を生成する．

ここで，Ỹ T は T 時点までのすべての観測値を表し，
β̃T = [β1 , β2 , . . . , βT ]γ̃T = [b11 , . . . , b1T , . . . , b41 , . . . , b4T , d11 , . . . , d1T , . . . , d41 , . . . , d4T ,

α11 , . . . , α1T , . . . , α41 , . . . , α4T ]θ̃T = β̃T , γ̃Tσ̃ 2 = [σβ2 , σb2 , σd2 , σα2 , σ 2 ]である．k は k = 1, 2, . . . , M1 , M1 + 1, . . . , M1 + M2 とし，M1 + M2 回の繰り返しの結果のうち，最初の M1 回の結果は捨て，残りの M2 回のサンプリングの結果から Θ の平均値を計算しそれを推定値とする．すなわち，E g(Θ|Ỹ T ) =

1
M2



M1 +M2g(Θ(k) |Ỹ T )k=M1 +1を推定値とする．ステップ
（a）
（
，b）
は具体的に以下の手順から構成される．
ステップ
（a）:
θt|t = E(θt |Ỹ t ), Σt|t = V (θt |Ỹ t ), θt+1|t = E(θt+1 |Ỹ t ), Σt+1|t = V (θt+1 |Ỹ t )とする．θ̃T は，



T −1

(4.4)

P (θ̃T |Ỹ T ) = P (θT |Ỹ T )P (θt |θt+1 , Ỹ t )t=1

によって生成する．(θT |Ỹ ), (θt |θt+1 , Ỹ ) は，
T

t

(4.5)

(θT |Ỹ T ) ∼ N (θT |T , ΣT |T )(4.6)

(θt |θt+1 , Ỹ t ) ∼ N (θt|T , Σt|T )

210

統計数理

第 51 巻

第 2 号 2003に従う．ここで，θT |T = E(θT |Ỹ T )ΣT |T = V (θT |Ỹ T )θt|T = E(θt |θt+1 , Ỹ T ) = E(θt |θt|t , θt+1 )Σt|T = V (θt |θt+1 , Ỹ T ) = V (θt |θt|t , θt+1 )である．(4.4) は θ̃T のサンプリングを，まず P (θT |Ỹ T ) から θT を生成し，その後 t = T − 1
から t = 1 まで逐次的に P (θt |θt+1 , Ỹ T ) から θt を生成することを示している．(4.5) における
θT |T ，ΣT |T は基本的にカルマン・フィルターのアルゴリズムにより得られる．与えられた初期値 θ0|0 = E(θ0 )，Σ0|0 = V (θ0 ) および σ̃02 に対し，以下の更新アルゴリズムにより逐次 θt|t ，Σt|tを計算する．
(i)

θt|t−1 = Ft θt−1|t−1
Σt|t−1 = Ft Σt−1|t−1 Ft + Qt ;(ii)

Qt = V (Jt ¯t ) = Jt V (¯t )Jt

jグループ j（j = 1, . . . , 8）の観測値 Ȳtj に対し，N (0, σ 2 )[vit≥ −Ht(j) θt|t−1 ] より，Nj Ȳtjj
個の vit を発生させる．ここで，Ht(j) は行列 Ht の第 j 行を表す．同様に N (0, σ 2 )[vit≤−Ht(j) θt|t−1 ] より Nj (1 − Ȳtj ) 個の乱数を発生させる．
(iii)

j∗
j
j∗
合計 Nj 個の vit に対し，yit= vit+ Ht(j) θt|t−1 とし，それらの yitの平均 Ȳtj∗ を求め，Ȳt∗ = [Ȳt1∗ , Ȳt2∗ , Ȳt3∗ , . . . , Ȳt8∗ ] とする．

(iv)

ηt|t−1 = y ∗it − Ht θt|t−1

ft|t−1 = Ht−1 Σt|t−1 Ht−1
+ R̄

θt|t = θt|t−1 + Kt ηt|t−1
Σt|t = Σt|t−1 − Kt Ht Σt|t−1−1とする．ここで，Kt = Σt|t−1 Ht ft|t−1である．

以上の手順に従い，θ1|1 , . . . , θT |T ，Σ1|1 , . . . , ΣT |T を得る．
(v)

続いて平滑化を行なう．
θt|T = θt|t + Σt|t Ft (Ft Σt|t Ft + Qt )−1 (θt+1|t − Ft θt|t )Σt|T = Σt|t − Σt|t Ft (Ft Σt|t Ft + Qt )−1 Ft Σt|tにより，θT |T , θT −1|T , . . . , θ1|T を得た後，(4.6）式より (θt |θt+1 ，Ỹ t ) を発生する．これが平滑化後の状態の推定値となる．
: σ̃ 2 の成分 σi2 は以下の事後分布から生成される．
ステップ
（ b）
事前分布 σi2 |θi ∼ IG( ν20 , δ20 ) とし，(a) で得られた推定値を用い，事後分布 σi2 |θi , Ỹ t ∼ IG( ν21 , δ21 ) を計算する．ここにν1 = ν0 + T

である．また識別可能性のため σ 2 = 1 とし，σβ2 ，σb2 ，σd2 ，σα2 に対しては，事後分布のパラメー
タをδ1 = δ0 +

(θ(i)t+1|T − F(i)t θt|T )2として発生させる．ここで θ(i) ，F(i) は σi2 に対応する適当な行を示す．

労働力調査とローテーション・サンプリング

5.

211

擬似データによる失業率の推定

労働力調査結果は毎月総務省統計局から公表されており，そこには就業状態別人口や年齢別
就業者数，失業率などの情報がクロス集計結果として与えられているが，ローテーションに対応した調査回数別の数値は与えられておらず，それゆえ，現実のデータを用いて上記モデルの有効性を確かめることはできない．実際，公表される失業率もローテーション構造を考慮した推定が行なわれてはいない．
そこで以下では完全失業率の時系列データをもとに，グループ別失業率の擬似データを作成
し，試験的に前節のモデルの適合を行なう．擬似データの作成方法は以下のようになされた．
1981 年から 2002 年まで（22 年間）の完全失業率の季節調整済み月次系列に対し，調査回数に応じて 0.6％，−0.4％，−0.2％，0.0％を平均とした −0.125％〜0.125％ の一様乱数を加える．バイアスの平均は 0 となっている．このようにして生成された調査回数別失業率の推移は図 3 のとおりである．8 グループの集計値データに対し，比較のため（1）ケース A：Bell and Carolan（1998）におけると同様，失業率の値 Ȳij そのものを複数の要因に分解することを考える．ここでは潜在変数は考えなくても良い．
（2）
ケース B：（4.1）〜（4.3）で記述されるモデルを用い変動要因を推定する，の 2 ケースの推定を行なった．ケース B においては N = 200 とし，N = 50 についても参考として分析を行なった．以上の結果は図 4〜図 6 に見られるとおりである．
ここでのモデルにおいてはいくつかの簡単化がなされている．たとえば，調査区効果は捨象されている．またローテーション・バイアスはランダム・ウォークではなく定数として扱われている．
図 4-a〜4-g はケース A における推定結果である．Gibbs Sampling の繰り返し回数は M1 =2000,
M2 =3000 とした．図 4-a はトレンド，4-b はバイアスの推定値を表す．トレンドはスムーズに変化している．バイアスの推定値は多少真の値と異なっている．また図 4-c〜4-f は M2 =3000
回の繰り返しにおけるパラメータや状態の推定値のヒストグラムを描いたものである．ヒスト
グラムはいずれも単峰であり，また σβ2 の推定値は定常な動きを示している．これらはいずれも繰り返しの結果，推定値が収束していることを反映している．図 4-g は失業率の予測値（平滑値）である．予測値は概ねデータの動きを追っている．

図 3.

調査回数別失業率．

212

統計数理

図 4-a.

図 4-b.

第 51 巻

第 2 号 2003

βt の推定値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

b1 , b2 , b3 , b4 の推定値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

図 4-c.

βt の推定値：t = 2002 年 12 月のヒストグラム．

労働力調査とローテーション・サンプリング

図 4-d.

図 4-e.

σβ2 の推定値：ヒストグラム．

σβ2 の推定値：iteration 回数でのプロット．

図 4-f.

2 の推定値：ヒストグラム．
σα

213

214

統計数理

図 4-g.

図 5-a.

第 51 巻

第 2 号 2003

失業率の予測値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

βt の推定値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

図 5-a〜5-g はケース B（標本数 N = 200）における推定結果である．トレンドの推定値はケー
ス A とやや異なるものの，よりスムーズに変化している．5-b のバイアスの推定値は 図 4-b と比べ変動が大きいが，バイアスの真の大きさの順序どおり推定されている．ただし，プロビッ
トモデルで推定される（4.2）の右辺は，個人の失業を規定する潜在変数の平均値であって集計値としての失業率ではない．それゆえ，バイアスの推定値は必ずしも現実の失業率に対して実験的に加えられたバイアスに対応しないことに注意したい．ケース A との違いは，図 5-g の失業
率の予測値
（平滑値）
である．推定値系列はケース A の結果と比べて，若干データの動きから離れている．これは，バイナリーなデータから得られる情報が必然的に小さいことを示唆している．ただし，図 5-c〜5-f の繰り返しにおけるパラメータや状態の推定値のヒストグラムは単峰であり，推定値が収束していることを示している．
また図 6-a〜6-d はケース B において，N =50 とした場合のトレンドの推定値，σβ2 の推定値
及び失業率の予測値である．プロビットによる定式化の場合，N =50 では Gibbs Sampling による推定結果が安定的ではないことがわかる．これは擬似データ数（N = 50）の問題であると同時にその質の問題でもある．データ数 N が小さいと，失業率が低いとき，擬似データの値

労働力調査とローテーション・サンプリング

図 5-b.

図 5-c.

b1 , b2 , b3 , b4 の推定値：時系列プロット．

βt の推定値：t = 2002 年 12 月のヒストグラム．

図 5-d.

σβ2 の推定値：ヒストグラム．

215

216

統計数理

図 5-e.

第 2 号 2003

σβ2 の推定値：iteration 回数でのプロット．

図 5-f.

図 5-g.

第 51 巻2 の推定値：ヒストグラム．
σα

失業率の予測値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

労働力調査とローテーション・サンプリング

図 6-a.

βt の推定値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

図 6-b.

b1 , b2 , b3 , b4 の推定値：時系列プロット．

図 6-c.

σβ2 の推定値：ヒストグラム．

217

218

統計数理

図 6-d.

図 7-a.

第 51 巻

第 2 号 2003

失業率の予測値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

βt の推定値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

図 7-b.

b1 , b2 , b3 , b4 の推定値：時系列プロット．

労働力調査とローテーション・サンプリング

図 7-c.

図 7-d.

219

σβ2 の推定値：ヒストグラム．

失業率の予測値：t = 1981 年 1 月〜2002 年 12 月の時系列プロット．

は多くの 0 を続けて発生することになる．すなわちほとんど生起していない現象を分析することになる．その結果推定値が大きく振れてしまう．しかし，ある地域における高齢者の失業率の推定などを行なうのは失業率が高いからであり，現実に 0 が続く事はないはずである．そこ
で，現実の失業率の値に 5％を加えることにより新たな擬似データを作成し，N = 50 として同様の実験を行なってみた．図 7-a〜7-d はその結果であるが，予測値には改善が見られよりデー
タの動きを追っている事がわかる．
以上のように，（4.1）〜（4.3）
で記述されるモデルはシミュレートされたデータに対しローテー
ションバイアスやトレンドなどを一定の精度で推定することに成功していることがわかる．小
さなデータ数の現実データに対する分析が待たれるところである．
6.

結語

本稿では，わが国の労働力調査を念頭におき，ローテーション構造をもつサンプリングにより得られた観測値を分析するための時系列モデルについて論じた．公表データにおいては，

220

統計数理

第 51 巻

第 2 号 2003

ローテーションに関する情報が得られないため，モデルを試験的に適合することに主眼がおかれた．もしわが国の労働力調査結果の集計表にローテーションに関する情報が附加されるならば，本稿で紹介したモデルの応用範囲は広いと思われる．以下にこのようなモデルの有用性について述べ結語としたい．
現在は標本の単純平均として失業率が計算されている．すでに述べたようにローテーション
構造を考慮していない推定は誤差分散が大きいはずである．本モデルを用いることにより，少なくともより有効な推定がなされうるはずである．
季節性要因の除去の問題も重要な実際的問題である．すでに述べたように，現実には曜日効
果やうるう年効果などを考慮に入れたより細やかな季節調整が好まれ，失業率の推定値を求めた後，X12-ARIMA 等によって季節調整値が計算されている．本モデルの状態変数として季節性要因を加えることにより，ローテーション構造が季節調整に与える影響を調べることは重要な課題であろう．
さらに本モデルを用いて調査方式の変更の影響を調べることも可能である．わが国の労働力調査の調査方式はたびたび変更されてきている．さらに IT 化の進行に伴い，今後インターネッ
トによる調査等も近々導入されるであろう．そのような場合に新たな調査方式のインパクトを把握することは重要な課題であろう．Bell（1998）は状態変数のなかにステップ関数を導入する
ことにより，電話によるインタビュー方式が労働力調査結果に与える影響を分析している．本
モデルを用いて，わが国における調査方式の変更の影響を調べることは是非試みるべき実験であろう．
しかし，モデルのより重要な応用は，ローテーション構造の決定であろう．わが国で用いられている 2-10-2（4），オーストラリアの労働力調査で用いられている 8-0-0（8）というローテーション構造は果たして妥当なものであろうか．どのようなローテーション構造の調査を行なう
べきか，またそもそもローテーション標本調査を行なうべきかは各国の現実のデータの振る舞
いに依存するはずである．上のモデルで言えば，σα2 や σb2 は標本調査にローテーションを導入するがゆえに生ずる誤差である．βt やその変化など，興味あるパラメータの変動の大きさとこ
れらの要因の相対的な大小関係によって，どのようなローテーション構造の調査を行なうべきか，そもそもローテーション標本調査を行なうべきか否か等が決定される．この点は，今後検討すべき重要課題である．
わが国の官庁統計を用いた実証研究はごく初歩的な分析例を除きあまり見られず，この点では諸外国に比べ大きく遅れている．それはわが国の官庁統計の一般利用が難しいことに起因す
る．今後はより多くの研究者に利用されるような調査結果の公表を考えるべきであろう．それによって研究者からのフィードバックがなされ，調査方法の改善も期待される．中長期的な視点に立ち，より質の高い統計の蓄積を目指すべきであろう．
謝辞

一橋大学経済学研究科の経済統計ワークショップ，東京大学経済学研究科の応用統計ワーク
ショップの参加者，および 2 人の匿名レフリーから貴重なコメントをいただいた．記して感謝致したい．なお，本研究は統計数理研究所共同研究プログラム（14-共研-2024），および文部科
学省科学研究費補助金（13553001, 14208023）の援助を受けている．

労働力調査とローテーション・サンプリング

221

参 考 文 献Bailar, B. A.（1975）. The eﬀect of rotation group bias on estimates from panel surveys, J. Amer.
Statist. Assoc., 70, 23–30.
. Trend estimation from small areas from a continuing survey with controlled sample
Bell, P. A.（1998）
overlap, Australian Bureau of Statistics, Working Paper, No. 98/1.
Bell, P. A. and Carolan, M.（1998）. Using state space models and composite estimation to measure
the eﬀects of telephone interviewing on labour force estimates, Australian Bureau of Statistics,
Working Paper, No. 98/2.
McLaren, C. H. and Steel, D. G.（1997）. The eﬀect of diﬀerent rotation patterns on the sampling
variance of seasonal and trend ﬁlters, Proceedings of the Survey Research Methods Section,
American Statistical Association, 790–795.
Kim, C-J. and Nelson C. R.（1999）. State-space Models with Regime Switching, The MIT Press,
Cambridge, Massachusetts.
Pfeﬀermann, D., Feder, M. and Signorelli, D.（1998）. Estimation of autocorrelations of survey errors
with application to trend estimation in small areas, J. Bus. Econom. Statist., 16, 339–348.
総務庁統計局（1999）. 『労働力調査標本設計の解説』．

222

Proceedings of the Institute of Statistical Mathematics Vol. 51, No. 2, 199–222 (2003)

Labor Force Survey and Rotation Sampling
Satoru Kanoh(Institute of Economic Research, Hitotsubashi University)

The labor force survey in Japan has a unique characteristic in its sampling scheme
called the rotation sampling. This paper proposes a time-series model for analyzing theresult of this survey taking this sampling scheme into consideration. In the time seriesmodel, individual unemployment is expressed by the probit model and its latent variable
is related to trend, group eﬀect and individual eﬀect. Because of the model’s non-linearity,
the maximum likelihood method using the Kalman Filtering algorithm is not applicablefor estimating the model’s parameters. In this paper, the estimation is carried out by the
Gibbs sampling technique. Unfortunately, the survey result does not include information
regarding the rotation groups and, for this reason, an experiment to validate the model
must rely on simulated data. Finally, a couple of important applications and extensionsof the model are given to conclude this paper.

Key words: Labor force survey, rotation sampling, probit model, Gibbs sampling, state spacepresentation.