ryu9_2007_2

ドキュメンテーションです。
プレゼンテーションでもいいです。
とりあえず，どっちでもいいです。

とりあえず与えられたテーマに沿って，何かコンテンツ（というほどのものでもないかもしれないが）を作って公開してもらいます。
これまでにやった解析手法を大いに活用してもらいたい。

１．テーマ設定
２．資料収集
３．データ作成
４．データ解析
５．結果の解釈
６．結果を『他人に理解してもらえるような形で記述』
７．公開

因みにWORD，Excelで作成したものは，そのままHTML形式で保存して，Web表示できます。
でかいです。
とてつもなくでかいですし，迷惑なことに，図形などは新たにディレクトリを作成してくれたりします。
ということはPowerPoint（以下PPT）でも同様です。

今回は，BlogとPPTを使ってコンテンツ（というほどのものでもないかもしれないが）を完成させます。
ところでこの教室はOffice2007ですが，私は持っていませんし，いまひとつ気が進みません。現在2003からOpenOffice.orgに移行中なのでどうしたものかと。
そこらへんは考えておきます。

△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼
とりあえず，習うより慣れろで，参考に一つ作ってみます。まねしてください
頭ではなく，手で理解してください。
「分かる」 ≠ 「できる」です。
△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼△▼
１．昨年の選挙で仲井間知事は失業率云々を公約にあげたが，失業率を考える。
２．とりあえず沖縄と全国の失業率を比べてみる。まずは沖縄県統計課へ。
３．失業率のデータを作ろう。ExcelをRで読み込めるように整形しよう。
４．時系列データなので時系列分析で将来予測をやってみるか。
５．なるほど，おお，未来が見えた？
６．とりあえずドキュメント作成。
７．公開

☐◆☐◆☐◆☐◆☐◆☐◆☐◆☐◆☐◆
ここから実際の作業をシミュレートしましょう
☐◆☐◆☐◆☐◆☐◆☐◆☐◆☐◆☐◆

沖縄県統計課に行く

　　　　　　　　　　　　　　　　　　　　　　　　　　どうも，これっぽい。
全国のデータは

　　　　　　　　　　　　　　　　　　　　　　　　　　どうも，それっぽい。
ここから「総務省」→「労働力調査」→「統計表一覧」→「長期時系列データ」をたどっていく。

両方のサイトからExcelのファイルを引っ張ってきて，整形する。

こういうときExcelはは不便だ。ASCII保存してAWKとかで編集できると楽。

メモ帳などのテキストエディタにExcelのデータを貼り付ける。
（もちろん，CSVで保存したり，直接読み見込むことも可能，例えば
こっちとあっち
特にあっちは秀逸

・*:.。..。.:*・∽・*:.。..。.:*・∽・*:.。..。.:*・∽・*:.。..。.:*・・*:.。..。.:*・∽・*
参考：　整形したExcel File 　　⇒　　からコピペしたテキストファイル
あくまで参考なので，自分で作業したほうが後々のため良い。
ただしExcel使いにだけはなりたくない，という志の高い方は
私が作ったやつをそのまま使ってください。
・*:.。..。.:*・∽・*:.。..。.:*・∽・*:.。..。.:*・∽・*:.。..。.:*・・*:.。..。.:*・∽・*
undat<-read.table("c:/00dat/unemp.txt",T)
plot(undat)
ts.plot(undat)


ts.plot(undat,gpars=list(xlab="年",ylab="失業率"),lty=c(1:2),col=c(1:2))
legend(locator(1),c("沖縄","全国"),lty=c(1:2),col=c(1:2))
##　凡例を置きたい場所をクリックする

多変量時系列分析をする訳ではないので，データ系列を分割しよう。
時系列オブジェクトに変換する前に
attchで扱ってもいいと思うが，とりあえず
変数名が日本語なので分割して英字にしよう。

attach(undat)
okn<-沖縄
jpn<-全国

okn<-ts(okn,start=c(1996,1),frequency=12)
jpn<-ts(jpn,start=c(1996,1),frequency=12)

## 勿論こっちでも良い。
## okn<-ts(undat$沖縄,start=c(1996,1),frequency=12)
## jpn<-ts(undat$全国,start=c(1996,1),frequency=12)

plot(okn)

横軸のラベルが多少良くなった。

時系列分析いきます。

（1） 自己相関，相互相関
par(mfrow=c(3,1))
acf(okn)
acf(jpn)
ccf(okn,jpn)

（2） 偏自己相関
par(mfrow=c(2,1))
pacf(okn)
pacf(jpn)

（3） ピリオドグラムとその平滑化
par(mfrow=c(2,2))
spec.pgram(okn)
spec.pgram(okn,spans=c(3,3))
spec.pgram(jpn)
spec.pgram(jpn,spans=c(3,3))

（4） スペクトル
par(mfrow=c(1,2))
spectrum(okn,method="ar")
spectrum(jpn,method="ar")

(3)はいらないと思うけど,(4)よりデータには
トレンドがあり非定常である。
なんのこっちゃ，
と思うかもしれませんが，スペクトルは周波数という視点から，データを眺めたもので，横軸は【１／ｆ】で周波数分の逆数＝周期を意味します。
ので，このデータには長周期あるいはトレンドが見出せるわけです。

スペクトルの見方は面倒なので，Web上にあげませんから。実習室でしっかり聞いといてください。
トレンド，周期成分，ホワイトノイズの３つを図示説明します。

（＿）/　とりあえず頭を空にして以下のコードを走らせてみよう
## １回差分をとってみる
par(mfrow=c(2,1))
acf(diff(okn))
acf(diff(jpn))

par(mfrow=c(2,1))
pacf(diff(okn))
pacf(diff(jpn))

par(mfrow=c(1,2))
spectrum(diff(okn),method="ar")
spectrum(diff(jpn),method="ar")

## ARIMA(1,1,1)にしてみる
par(mfrow=c(1,1))
(okar<-arima(okn,order=c(1,1,1)))
ts.plot(okar$residuals)
(okpr<-predict(okar,n.ahead=24))
okup<-okpr$pred+1.96*okpr$se
okdn<-okpr$pred-1.96*okpr$se
ts.plot(okn,okpr$pred,okup,okdn)

(jpar<-arima(jpn,order=c(1,1,1)))
ts.plot(jpar$residuals)
(jppr<-predict(jpar,n.ahead=24))
jpup<-jppr$pred+1.96*jppr$se
jpdn<-jppr$pred-1.96*jppr$se
ts.plot(jpn,jppr$pred,jpup,jpdn)

par(mfrow=c(2,1))
ts.plot(okn,okpr$pred,okup,okdn,gpars=list(lt=c(1,2,3,4),col=c(1,2,3,4)))
ts.plot(jpn,jppr$pred,jpup,jpdn,gpars=list(lt=c(1,2,3,4),col=c(1,2,3,4)))

##　安易な予測を反省してみる
##　ARIMA(10,1,2)！！で盲省する
(okar<-arima(okn,order=c(10,1,2)))
(okpr<-predict(okar,n.ahead=24))
okup<-okpr$pred+1.96*okpr$se
okdn<-okpr$pred-1.96*okpr$se

##　ARIMA(10,1,10)！！で脱力する
(jpar<-arima(jpn,order=c(10,1,10)))
(jppr<-predict(jpar,n.ahead=24))
jpup<-jppr$pred+1.96*jppr$se
jpdn<-jppr$pred-1.96*jppr$se

par(mfrow=c(2,1))
ts.plot(okn,okpr$pred,okup,okdn,gpars=list(lt=c(1,2,3,4),col=c(1,2,3,4)))
ts.plot(jpn,jppr$pred,jpup,jpdn,gpars=list(lt=c(1,2,3,4),col=c(1,2,3,4)))


SAMPLEのBlogはここ，PPTはここにあります。

本日はとりあえずここまでとし，続きは次回とします。

皆さんは，統計データのサイトやらをブラウジングして
扱ってみたいデータとかを探しておいてください。