時系列分析と状態空間モデルの基礎：サポートページ
『時系列分析と状態空間モデルの基礎：RとStanで学ぶ理論と実装』のサポートページです。
本書に使用したサンプルデータとR,StanのコードはすべてGitHubから参照できます。
緑色の「Clone or download」というボタンをクリックしてから「Download ZIP」をクリックすると、すべてのファイルをZIP形式でダウンロードできます。
書籍のサンプルコードとデータ
注意事項（2019年7月14日追記）
Nipponパッケージが使えなくなってしまったため、p273の「is.jholiday」関数が使えなくなってしまいました。
同様の機能を持つ関数を当方が作成しました。Nipponパッケージを読み込む『library(Nippon)』の代わりに、以下のコードを実行してください。祝日判定が可能になります。
詳細は『R言語における日本の祝日判定』を参照してください。
# 関数の読み込み
source("https://raw.githubusercontent.com/logics-of-blue/website/master/010_forecast/20190714_R%E8%A8%80%E8%AA%9E%E3%81%AB%E3%81%8A%E3%81%91%E3%82%8B%E6%97%A5%E6%9C%AC%E3%81%AE%E7%A5%9D%E6%97%A5%E5%88%A4%E5%AE%9A/jholiday.R", encoding="utf-8")
注意事項（2019年1月21日追記、1月26日再修正）
新しい『forecast_8.4』のバージョンにおいては、auto.arima関数の実行結果が、本書第2部第7章p107と異なるものとなります。これはforecastパッケージのバージョンが上がり、モデル選択のデフォルト設定が変わったことが理由です。p106～p107にかけてのauto.arima関数の実装コードに『seasonal.test = “ch”』または『seasonal.test = “ocsb”』という指定を追加することで、本書と同じ結果が得られます。修正後のコードは以下のようになります。
sarimax_petro_law <- auto.arima( y = train[, "front"], xreg = petro_law, ic = "aic", seasonal.test = "ch", max.order = 7, stepwise = F, approximation = F, parallel = T, num.cores = 4 )
(2019年5月6日追記)
新しい『forecast_8.7』のバージョンにおいては、p112の予測のコードがエラーとなってしまいます。これはforecastパッケージのバージョンが上がり、外生変数の入れ方が変わってしまったことが理由です。
p112の実行コードを以下に変更すると、本書と同様の結果が得られます(as.matrix関数を適用しました)。
# 過去の石油価格の平均値を使う petro_law_mean <- data.frame( PetrolPrice=rep(mean(train[, "PetrolPrice"]),12), law=rep(1, 12) ) sarimax_f_mean <- forecast(sarimax_petro_law, xreg=as.matrix(petro_law_mean)) # 直近の石油価格を使う petro_law_tail <- data.frame( PetrolPrice=rep(tail(train[, "PetrolPrice"], n=1),12), law=rep(1, 12) ) sarimax_f_tail <- forecast(sarimax_petro_law, xreg=as.matrix(petro_law_tail))
(2019年2月26日追記)
新しい『prais 1.1.0』においては、prais.winsten関数がなくなりprais_winsten関数に名称が変更されています。引数にも一部変更があります。p136におけるパッケージを使ったPrais-Winsten法の実行コードを以下に変更すると、本書とほぼ同様の結果が得られます。
mod_gls_PW <- prais_winsten(y_ar ~ x_ar, data=d, max_iter = 1) summary(mod_gls_PW)
（2018年4月15日追記）
第5部11章(p272)のautoplot関数の実行時にエラーとなった場合は、RStudioを再起動したうえで、以下のコードを実行してみてください。
zooパッケージの仕様変更に伴い、グラフが正しく描画できなくなるエラーはこれで解決できるはずです。
install.packages("devtools")
require(devtools)
install_version("zoo", version = "1.8-0", repos = "http://cran.us.r-project.org")
（2019年1月21日追記）
第5部11章(p272)のautoplot関数の実行時にエラーに関する追加情報です。
以下のパッケージのバージョンにおいては、正しく動作することが確認されました。 パッケージのバージョンを最新にすることで、エラーは出なくなるはずです。
ggfortify_0.4.5・ggplot2_3.1.0・xts_0.11-2・zoo_1.8-4
（2018年10月1日追記）
出版社サイトにて、正誤表が公開されています。合わせて参照してください。
時系列分析と状態空間モデルの基礎： RとStanで学ぶ理論と実装  ハヤブサの表紙が目印です。 Twitterでは「#隼時系列本」のハッシュタグをお使いください。 Amazonでは在庫切れで価格が高騰することがあります。 |
 |
丸善さん/ジュンク堂書店さんの在庫はこちらから確認できます。
紀伊国屋書店さんの在庫はこちらから確認できます。
書店に置かれていなくても、取り寄せてもらうことができるはずです。
出版社のWebサイトから直接購入することもできます。
送料は無料で後払いです。
出版社のWebサイトはこちらです
基本情報
- 出版社 ： プレアデス出版
- 著者 ： 馬場真哉（このサイト、Logics of Blueの管理人です）
- タイトル ： 時系列分析と状態空間モデルの基礎：RとStanで学ぶ理論と実装
2018年2月中旬発売。
2018年4月に重版となりました(2018年4月15日追記)。
2018年10月に再度重版となりました。3刷目となります(2018年10月1日追記)。
2019年3月に再度重版となりました。4刷目となります(2019年5月6日追記)。
かなり詳細に目次やR関数一覧を作っています。
書籍を購入されるかどうか、こちらを見てご判断いただいても良いかと思います。
2018年2月7日追記
本文の1部を公開しました。状態空間モデルを推定するためのパッケージであるKFASの使い方に関する内容です。
・第5部8章ー実装：KFASの使い方
2018年4月15日追記
出版社のサイトから、第1部の1~2章を試読することができるようになりました。
お勧め書籍｜プレアデス出版
書籍の特徴
Box-Jenkins法と状態空間モデルという、時系列分析の二大フレームワークを一気に学べる、時系列分析の入門書です。状態空間モデルは、線形ガウス状態空間モデルだけでなく、一般化状態空間モデルについても言及しました。
また、見せかけの回帰・VARモデル・GARCHモデルといった重要なトピックも併せて解説。
時系列分析を初めて学ぶという方でも、基礎から応用まで順を追って学べます。
全体を通して「時系列モデルの”考え方”」がつかめるように配慮しました。
そのため、第1部では分析ツールの話は出てこず、データ生成過程の解説から始まります。この方が、最終的には実務においても役に立つと信じています。
もちろんタイトルの通り、RとStanを用いた実装の仕方も解説しています。forecastやKFASパッケージを用いた効率的な分析、ggplot、ggfortifyによる美麗なグラフ、そしてStanを用いたベイズモデリングを学びます。
書籍の内容
タイトル「時系列分析と状態空間モデル」の通り、全体の半分が状態空間モデルで、残りの半分がBox-Jenkins法などの時系列モデルの解説となっています。
状態空間モデルの前座としてARIMAの紹介を挟んでいるわけではありません。
この本では、データの変換やARIMA次数の決定方法、季節性・外生変数の組み込み方、モデルの残差チェックなど、Box-Jenkins法を使いこなすための様々な技術を解説しています。
Box-Jenkins法は「時系列データをどのように分析するか」という知恵と工夫の集積です。Box-Jenkins法を学ぶことで、データの定常性や和分過程といった基礎概念や残差の自己相関のチェックといった実務的な技術を体系的に学ぶことができます。
もちろん状態空間モデルも指数を割いて解説しています。
状態空間モデルはやや複雑な概念です。そのため、状態空間モデルの導入・線形ガウス状態空間モデル・一般化状態空間モデルと3つの部に分けて解説を行いました。
線形ガウス状態空間モデルを学ぶパートでは、ローカルレベルモデルを対象として、カルマンフィルタと散漫カルマンフィルタ、対数尤度と散漫対数尤度を、ライブラリを用いずに実装します。dlmパッケージではカルマンフィルタが、KFASパッケージでは散漫カルマンフィルタが実装されています。両パッケージでなぜ推定結果が異なるのか、その理由がわかるはずです。
また、ローカル線形トレンドモデルなどは、推定された「トレンド」の解釈で戸惑う方も散見されます。1次のトレンドや2次のトレンドという用語を学ぶだけでなく、トレンドが変化するとはどのような状況か、そもそもトレンドとは何者か、数式・日本語・シミュレーションを通して丁寧に解説します。
応用編において用いるパッケージはKFASで統一しました。dlmに関しては当サイトで様々な解説がありますし、KFASの方が実行速度が速いことが理由です。1つのパッケージに絞ることで、覚えるべき内容が減り、状態空間モデルの”考え方”に注力できるようになります。非ガウシアンなモデルはStanに譲ることで、覚えるべき文法をさらに減らしました。
代わりに線形ガウス状態空間モデル推定におけるKFASの使い方については妥協なく解説しました。パラメタ推定・フィルタリング・平滑化・予測・補間・信頼区間と予測区間の使い分け・トレンドや季節成分の抽出など、この本を読めば、KFASの基本的な使い方に困ることはないはずです。
一般化状態空間モデルを学ぶパートでは、「ベイズの定理とは何か」という基本から解説をはじめます。ベイズ推論においてなぜ乱数生成アルゴリズムが必要なのか、そしてハミルトニアンモンテカルロ法はどのようにして効率的なサンプリングを可能としているのか、その概要をつかむことができます。
応用編ではランダムエフェクトの入ったポアソン分布に従うデータなどを扱います。Stanのインストールから一般化状態空間モデルの推定までしっかりと学べます。
大きなフレームワークからは外れますが、見せかけの回帰・VARモデル・GARCHモデルといった「時系列分析をするにあたって、是非知っておいてほしい考え方」については章を分けて解説しています。
「これはやっちゃダメ」とか「これを使うべき」とかいったルールを覚えるのではなく、そもそも「時系列データを分析するとは一体何か」を学んでほしい。
こういったことを念頭に、初めて時系列分析を学ぶ人に知っておいてほしい”考え方”を解説しました。
この本に載っていないこと
この本だけでは学ぶことができない内容もあります。
”考え方”の説明に注力した分、実際にスクラッチから時系列モデルを構築するような技術的な側面は大きく削りました。GARCHなどのパラメタ推定の方法についてはほぼノータッチです。
カルマンフィルタや平滑化についてもローカルレベルモデルしか解説していません。より一般的な線形ガウス状態空間モデルのカルマンフィルタの計算式などは成書に譲ることとなりました。
ベイズ推論についても同様で、あくまでもそのアイデアを述べるにとどまっています。
この本でもたびたび引用させていただいた、沖本(2010)『計量時系列分析』や野村(2016)『カルマンフィルタ』、松浦(2016)『RとStanでベイズ統計モデリング』、岩波データサイエンス刊行委員会(2017)『岩波データサイエンスVol6』といったより高度な書籍への橋渡しを狙って執筆したこともあります。隼時系列本を読んだ後にこれらの書籍を読まれると、より理解が深まるかと思います。
状態空間モデルにおいては、カルマンフィルタとMCMCによる推定方法しか記載がありません。
粒子フィルタなどの手法については言及しておりません。
また、以下の内容についても記載がありません
・facebook社のprophet
・機械学習法による時系列分析
・指数平滑化法やHolt-Winters法、X11といった統計モデル以前の潮流
データ生成過程を推定するという大きな思想から外れてしまい、全体のバランスをとることが難しいと考えたからです。
この本を読む前提となる知識
時系列分析は統計学の応用編です。
隼時系列本は、統計学に関する深い知識を要求しないように配慮しましたが、検定や推定、回帰分析や最小二乗法などはある程度知っておくと理解が深まります。
統計学の基礎に自信がないという方は、前著『平均・分散から始める一般化線形モデル入門』を読まれると良いかと思います。
『平均・分散から始める一般化線形モデル入門』の次に読む本として隼時系列本を提案したという経緯もありますので、難易度的にはちょうどいいはずです。表紙もほとんど同じデザインで作っていただいております。
もちろん内容に直接の関係はないですので、隼時系列本だけでも完結しています。
『平均・分散から始める一般化線形モデル入門』のサポートページ
更新履歴
2018年1月19日：新規追加
2018年2月07日：書籍本文を一部公開したため、その旨を追記。
2018年2月27日：書籍の購入方法について追記。
2018年3月09日：書籍のリンク先を変更
2018年4月15日：autoplot関数のエラー対応・出版社の試読ページへのリンク・重版した旨を追記
2018年4月25日：書籍のコードとデータをGItHubで管理するようにしたので、ダウンロード先のURLを修正
2018年10月1日：重版した旨と正誤表の情報を追記
2019年1月21日：パッケージのバージョンが上がったことに伴う注意事項を追記
2019年1月26日：季節差分の判定のための検定手法を追記
2019年2月26日：praisパッケージのバージョンが上がったことに伴う注意事項を追記
2019年5月06日：パッケージのバージョンが上がったことに伴う注意事項を追記・重版した旨を追記
2019年7月14日：Nipponパッケージが使えなくなったので、代替手段を追記
“時系列分析と状態空間モデルの基礎：サポートページ” に対して114件のコメントがあります。
コメントを残す
このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください。
新刊が出るのですね! 早速アマゾンで予約しました。楽しみです。
ありがとうございます！
購入を検討しているのでが都合により電子書籍形式での販売を希望しています。
電子書籍化の予定はありますか？
申し訳ないのですが、当面の間は電子版は出ないものと思います。
できればそのうち出したいなと思うのですが。。。
素晴らしい書籍ありがとうございます。購入して気になる章から読み解いております。早速ですがp.272のautoplotの行を実行すると、下記エラーが帰ってくるのですが何か要因ございますでしょうか？
R初心者でつまらない質問で申し訳ございませんが、お手すきの時にでもご確認頂ければ幸いです。よろしくお願いします。
Error in data.frame(index(model), …) :
arguments imply differing number of rows: 77, 1, 0
書籍をお読みいただきありがとうございます。
管理人の馬場です。
お問い合わせの件ですが、私の方で実行すると、エラーは出ませんでした。
おそらく実行環境に問題があるのではないかと思います。
エラーメッセージを見ると、データフレームの長さが合わないと怒られているので、別のデータが使われているような気がします。
RStudioをいったん落としてから再度実行すればうまくいくかもしれません。
ところで、Rのバージョンは3.4.3をお使いでしょうか。
また、念のため以下のコードを実行してライブラリを最新にしてください。
install.packages(“KFAS”)
install.packages(“xts”)
install.packages(“Nippon”)
install.packages(“ggplot2”)
install.packages(“ggfortify”)
install.packages(“gridExtra”)
そのうえで、RStudioをいったん右上の×ボタンを押して消した後、再度起動してから、書籍本文のコードを実行してください。
あるいは、この記事の上部で配布している『tsa-ssm-book-data.zip』の『5-11-応用：周期性のある日単位データの分析.R』を実行してください。
ご確認お願い致します。
読みやすい本をありがとうございます。とても細かいのですが、38ページの(2-6)式に「ホワイトノイズ＋定数も定常過程であるといえます。」とありますが、(2-6)式は「ホワイトノイズ＋1個ずらしたホワイトノイズ＋定数」であって「ホワイトノイズ＋定数」ではないのではないでしょうか。1次の自己相関がゼロにならないと思います。いずれにせよ定常過程ではあるのですが。（と思ったのですが、勘違いであったら申し訳ありません。）
C. K.様
書籍を丁寧にお読みいただき、ありがとうございます。
1階差分でトレンドを消せる、という説明の箇所ですね。自己相関が残るという点はその通りです。
説明の意図としては「定常過程になった」ということだけに着目しているので、記述として誤りではないと思います。
ただ、おっしゃる通り「ホワイトノイズ＋1個ずらしたホワイトノイズ＋定数」と書くのが丁寧ではありますね。ホワイトノイズの和は定常過程なので、1個ずらしたというところには着目していませんでした。
あと、違和感を覚える理由として、式(2-5)のようなデータには普通、差分はとらないというのもありますね。時間を説明変数とした線形回帰モデルで表現できますので。
あまり良い例ではなかったかもしれません。
(誤りとまではいかないという認識ですので)修正はかけないかと思いますが、ご指摘ありがとうございます。
もし私の認識違いがあれば、遠慮なくご指摘ください。
早速のご返信ありがとうございます。ご説明ですっきりしました。本筋ではない細かいご質問に丁寧に回答いただきありがとうございました。
2018年2月27日 at 10:47 AMの匿名希望さんと同様のエラーが発生してしまいます。
フィルタリング、スムージングまで実施できるのですが、
salesのautoplotだけ実行できません。
再起動しても、ご指定のファイルのコードでもできませんでした。
ご確認お願いします。
YA様
管理人の馬場です。
コメントありがとうございます。
ご不便おかけしております。
こちらで実行したら、やはり動きます。
再起動というのは、RStudioをいったん落とし、再度立ち上げてから、その直後に、
第5部11章のコードのみを実行したということで間違いないでしょうか。
また、他の章でのautoplotは動きますでしょうか。
autoplotはすべてほぼ同じコーディングとなっているので、
5部11章だけ動かないというのは、不思議です。
こちらも再起動の上で実行してみて、エラーが出るかどうか確認してみてください。
また、autoplot以前のコード、例えば
file_data <- read.csv("5-11-sales_data.csv")
sales <- as.xts(read.zoo(file_data))
の段階では、エラーは出ていないでしょうか。
お二人から同様のエラーが出た以上、何らかの原因があると思うのですが、
こちらではなかなか再現できません。
お手数ですが、YA様の実行環境とエラーが出た再現手順を教えてもらえますでしょうか。
なお、当方の動作環境は、p72に記載のある通りです。
よろしくお願いします。
ご連絡が遅くなり失礼しました。
再起動および第5部11章のコードのみ実行すると、
以下の三箇所のみエラーが発生します。
> autoplot(sales, main = “架空の売上データ”)
Error in data.frame(index(model), …) :
arguments imply differing number of rows: 77, 1, 0
> p_data grid.arrange(p_data, p_trend, p_cycle)
Error in arrangeGrob(…) : object ‘p_data’ not found
その他は、他のautoplot含めすべて実行できています。
お手数をおかけしますがご確認お願いいたします。
すみません、
もう一箇所は以下です。
> p_data <- autoplot(sales, main = "元データ")
Error in data.frame(index(model), …) :
arguments imply differing number of rows: 77, 1, 0
YA様
ありがとうございます。
以下の点、追加で教えていただけますでしょうか。
①
『autoplot(sales, main = “架空の売上データ”)』でエラーが出るまでの、実行されたコードとその結果を、最初からすべて載せてもらえますか？
コンソールをコピペしてもらえればと思います。
②
以下のコードの実行結果を教えてもらえますか？
コンソールをコピペしてもらえればと思います。
sessionInfo()
お手数ですが、お願いいたします。
①
【コード】
library(KFAS)
library(xts)
library(Nippon)
library(ggplot2)
library(ggfortify)
library(gridExtra)
setwd(“C:\\Users\\******”)
file_data <- read.csv("5-11-sales_data.csv")
head(file_data, n = 3)
sales library(KFAS)
> library(xts)
要求されたパッケージ zoo をロード中です
次のパッケージを付け加えます: ‘zoo’
以下のオブジェクトは ‘package:base’ からマスクされています:
as.Date, as.Date.numeric
> library(Nippon)
要求されたパッケージ maptools をロード中です
要求されたパッケージ sp をロード中です
Checking rgeos availability: FALSE
Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
which has a restricted licence. It is disabled by default;
to enable gpclib, type gpclibPermit()
要求されたパッケージ stringr をロード中です
要求されたパッケージ sf をロード中です
Linking to GEOS 3.6.1, GDAL 2.2.0, proj.4 4.9.3
> library(ggplot2)
> library(ggfortify)
> library(gridExtra)
>
> setwd(“C:\\Users\\******”)
> file_data head(file_data, n = 3)
date sales
1 2010-03-01 10
2 2010-03-02 34
3 2010-03-03 18
>
> sales head(sales, n = 3)
[,1]
2010-03-01 10
2010-03-02 34
2010-03-03 18
>
> autoplot(sales, main = “架空の売上データ”)
Error in data.frame(index(model), …) :
arguments imply differing number of rows: 77, 1, 0
②
> sessionInfo()
R version 3.4.3 (2017-11-30)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1
Matrix products: default
locale:
[1] LC_COLLATE=Japanese_Japan.932 LC_CTYPE=Japanese_Japan.932 LC_MONETARY=Japanese_Japan.932
[4] LC_NUMERIC=C LC_TIME=Japanese_Japan.932
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] gridExtra_2.3 ggfortify_0.4.2 ggplot2_2.2.1 Nippon_0.6.5 sf_0.6-0 stringr_1.3.0
[7] maptools_0.9-2 sp_1.2-7 xts_0.10-1 zoo_1.8-1 KFAS_1.3.1
loaded via a namespace (and not attached):
[1] Rcpp_0.12.15 urca_1.3-0 pillar_1.1.0 compiler_3.4.3 plyr_1.8.4 bindr_0.1
[7] class_7.3-14 vars_1.5-2 tools_3.4.3 tibble_1.4.2 nlme_3.1-131 gtable_0.2.0
[13] lattice_0.20-35 pkgconfig_2.0.1 rlang_0.1.6 DBI_0.8 bindrcpp_0.2 e1071_1.6-8
[19] dplyr_0.7.4 classInt_0.1-24 lmtest_0.9-35 grid_3.4.3 glue_1.2.0 R6_2.2.2
[25] foreign_0.8-69 udunits2_0.13 tidyr_0.8.0 purrr_0.2.4 magrittr_1.5 units_0.5-1
[31] scales_0.5.0 assertthat_0.2.0 colorspace_1.3-2 sandwich_2.4-0 stringi_1.1.6 lazyeval_0.2.1
[37] munsell_0.4.3
ディレクトリのパスはマスキングさせていただいています。
ご確認お願いいたします。
head(sales, n = 3)
autoplot(sales, main = “架空の売上データ”)
【実行結果】
> library(KFAS)
> library(xts)
なぜか12行目と13行目の間が省略されてしまっていたので、上記追記お願いします。
よろしくお願いします。
YA様
当方でもエラーを再現できました。
問題があったのはzooパッケージです。
当方の動作確認の際、xtsを最新にしたんですが、zooは最新になっていませんでした。
zooを最新のバージョンにすると、このエラーが出るようです。
というわけで、いったんRStudioを再起動してもらってから、以下のコードを実行してください。
install.packages(“devtools”)
require(devtools)
install_version(“zoo”, version = “1.8-0”, repos = “http://cran.us.r-project.org”)
これで、当方の執筆時のzooのバージョンとなります。
この後でコードを実行すれば動くはずです。
ご確認をお願いいたします。
同じエラー、現象で悩んでいたのですがzooのバージョンを指定してインストールしたら解消しました。
ありがとうございました。
馬場様
この本を購入し、まずはRによる時系列の勉強が進むと楽しみです。
さて、p.90に記載のある、”5-2-1-timeSeries.csv”のありかがサポートページを見てもよくわかりませんでした。
お手数ですが、ご教示をお願いいたします。
T.S様
ご購入いただき、ありがとうございます。
サンプルデータですが、サポートページ上部の
『サンプルデータとコードをダウンロード tsa-ssm-book-data.zip』
というリンクをクリックしていただければ、ダウンロードできるかと思います。
お手数ですが、ご確認をお願いいたします。
早速のご回答ありがとうございます。
これで先に進むことができます。
ありがとうございました。
馬場様
この度、本を購入させて頂きました。
以前の一般化線形モデルから大変分かり易く、数多ある統計解析の本でも記憶に残る素晴らしい本であると思っております。
一つご質問になります。
P106のauto.arima関数ですが、書籍の通り同様に計算を行っておりますが、P108のような結果となりませんでした。
Series: train[, “front”]
Regression with ARIMA(2,0,1)(0,1,1)[12] errors
Coefficients:
ar1 ar2 ma1 sma1 PetrolPrice law
1.1225 -0.1322 -0.8690 -0.8183 -0.3748 -0.3431
s.e. 0.0906 0.0876 0.0443 0.1129 0.1000 0.0473
sigma^2 estimated as 0.007624: log likelihood=168.12
AIC=-322.23 AICc=-321.53 BIC=-300.36
このように何かが原因で次数は計算の都度で異なってくるものなのでしょうか？
もし可能でしたら何かアドバイスなどご教示頂ければ幸いと存じます。
宜しくお願い致します。
R.Y様
書籍をお読みいただきありがとうございます。
次数ですが、forecastパッケージの2018年4月リリース版(V8.3)のバグかと思います(書籍では8.2を使っていました)。
parallel=F,trace=Tにして確認すると、差分の次数が正しく得られていないようでした。
現在のバージョンで正しく次数を選ぶためには『d=1,D=0』を引数に追記する必要があります。
こちらの方がAICは小さくなります。
あるいは、以下のコードを実行して、1つ前のバージョンに戻す方法もあります。
install.packages(“devtools”)
require(devtools)
install_version(“forecast”, version = “8.2”, repos = “http://cran.us.r-project.org”)
お手数ですが、ご確認をお願いいたします。
馬場様
お世話になります、R.Yです。
頂いた内容で結果の一致を確認することができました、ありがとうございます。
『d=1,D=0』とのことで階差をあえて指定しているといったことでしょうか。
当面はバージョンを戻すことで、対応してゆきたいと思っております。
お忙しいのにも関わらず、サポート頂きまして、ありがとうございました。
7月の大井町でのセミナー、お会いできること心より楽しみにしております。
宜しくお願い致します。
R.Y様
うまく動いたようでよかったです。
こちらこそ、どうぞよろしくお願いいたします。
p36の標本自己共分散は、k個サンプルをシフトして共分散を計算していくので、
1/nΣ…
ではなく
1/(N-k)Σ…
が正しいでしょうか？
hikosan様
コメントありがとうございます。
難しいところですが、ここはnをそのまま使うことが多いですね。
沖本(2010)など多くの書籍でも、同様となっています。
こちらの方が計算結果が安定して使いやすいということが理由です。
ありがとうございます。沖本本に注記がありましたね。ありがとうございました。
間違いだったらすみません。
P49の2-24式の次の行の
y_t-2も同様にφ_1*y_t-3+ε_t-3で表現できる
→ε_t-2の誤植？
2-25式
Σi=1 to m-1が正しい？
hikosan様
コメントありがとうございます。
ご指摘の通り、これは誤植かと思われます。
まずは書籍にミスがあったこと、失礼いたしました。
また、ご指摘をいただけましたこと、感謝いたします。
再度内容確認ののち、出版社と対応を検討いたします。
ありがとうございます。
ご確認ありがとうございました。
馬場様
とても読み易く役立つ本であり、とても感謝しております。
ところで、P.289の条件付確率の説明ですが、P(θ=女性|y=赤鞄発見)なので、分母は「女性である」10ではなく、「赤鞄所持」の4でしょうか。次の同時確率は、P(y=赤鞄発見|θ=女性)を使って説明してあるため、条件付確率の式が誤記なのかもしれませんが。
mtakayuki様
管理人の馬場です。
コメントありがとうございます。
また、拙著をお読みいただきありがとうございます。
該当箇所ですが、確かに矛盾した表現が見られますので、
誤植の修正も含めて対応を検討します。
まずは書籍にミスがあったこと、失礼いたしました。
また、ご指摘をいただけましたこと、感謝いたします。
馬場様
この度,本を購入させていただきました。
時系列を自分なりに学び難航しているなか,この本には何度も助けていただきました。ありがとうございます。
質問できる人がいないためここにコメントさせていただくのですが,ARIMAモデルでは異常値を表すことが出来ず,ARIMAXモデルはその欠点を外生変数で補っているという解釈で間違いありませんでしょうか？
お手数ですが解答のほどよろしくお願い致します
K様
コメントありがとうございます。
管理人の馬場です。
書籍を購入いただきありがとうございます。
お役に立てたようであれば幸いです。
質問に関してですが、「異常値」が指す内容が様々あるので、
一言では答えにくいです。
ARIMAでは外部のイベントに対応できず、
ARIMAXならそれに対応できる、
という意味ならば、正しい解釈だと思います。
異常値が何らかの外部のイベントによって発生する場合は、
ARIMAXで対応できるかもしれません。
理解することができました。
迅速なご返答ありがとうございました。
馬場様
KFASでSSMregressionにより外生変数を取り込んだ際の、
予測がうまくいきません。newdataの置き方が間違っていると考えているのですが、
アドバイス頂けないでしょうか。
恐れ入りますがよろしくおねがいします。
library(KFAS)
library(xts)
library(Nippon)
library(ggplot2)
library(ggfortify)
library(gridExtra)
library(Metrics)
nof<-read.csv("sale.csv",header=TRUE)
nof<-as.xts(read.zoo(nof))
nof_train<-nof[1:800,]
nof_test<-nof[801:805,]
aa<-nof_train[,2]
bb<-nof_train[,3]
build_kfas<-SSModel(H=NA,as.numeric(nof_train[,1]) ~SSMtrend(degree=1,Q=NA)+SSMseasonal(period=12,sea.type="dummy",Q=NA)+SSMregression(~as.numeric(aa)+factor(bb),Q=NA,type="common"),data=nof)
fit_kfas<-fitSSM(build_kfas,inits=c(1,1,1,1))
result_kfas<-KFS(fit_kfas$model,filtering=c("state","mean"),smoothing=c("state","mean"))
mu_filter_kfas<-result_kfas$a[-1]
mu_smooth_kfas<-result_kfas$alphahat
aa1<-nof_test[,2]
bb1<-nof_test[,3]
newdata1<-SSModel(formula=rep(NA,5)~SSMtrend(degree=1,Q=NA)+SSMseasonal(period=12,sea.type="dummy",Q=NA)+SSMregression(~as.numeric(aa1)+factor(bb1),Q=NA,type="common"),data=nof)
forecast_pred<-predict(fit_kfas$model,newdata=newdata1,type="response",interval="prediction",level=0.95)
Ichi様
コメントありがとうございます。
管理人の馬場です。
コードを見る限り、書籍に関する内容ではなく、純粋にKFASについての質問ですね。
データがないので推測になりますが、コードを見て気になったところを回答します。
1．パラメタが取得されていない
下記のように、最尤法で得られたパラメタを取得する必要があります。
pars <- exp(fit_kfas$optim.out$par) そのうえで、以下のように、各分散の値をparから引っ張ってくる必要があります。
newdata1<-SSModel( H = pars[4]。。。以下略 推定された分散は4つあるので、各々parsの中身を入れておきましょう。
2．応答変数の数がおかしい
NAの数が5つしかないのは間違いだと思います。
以下のように、テストデータの長さと一致させる必要があります。
newdata1<-SSModel( H = pars[4], formula= rep(NA,nrow(nof_test))。。。以下略
馬場様
返信頂きまして、誠にありがとうございます。
書籍の第5部-11章を基に、実践させて頂いています。
ご指摘いただいた箇所を修正したのですが、まだ下記のエラーが出ます。
Error in is.SSModel(newdata, na.check = TRUE, return.logical = FALSE) :
System matrices (excluding Z) contain NA or infinite values,
covariance matrices contain values larger than 1e+07
データの説明を簡単にさせて頂くと、目的変数のところが販売量、
説明変数aaが最高気温、bbが曜日となっています。また、読み込むcsvファイルの1列目を日付としています。
度々で申し訳ございませんが、よろしくお願いします。
馬場様
大変申し訳ないのですが、誤って本名を書いてしまいました。
名前のところだけ、変更して頂けますか。
お手数をお掛けします
対応しました。
馬場様
すみません、ご指摘頂いたnewdata1のH=parsの中身の入れ方が
よくわからないので、恐れ入りますが、説明して頂けないでしょうか。
（H=pars[1:4]のようにしてもうまくいきません。）
お手数ですが、よろしくお願いします。
H = pars[4]
のように、1つだけを入れます。
ちなみにparsの最後の1つが観測誤差（Hに入る）です。
parsが4つしかなかったら『H = pars[4]』ですね。
Qもparsの値を入れてください。
pars[1]とか pars[2]とか言ったように入れていきます。
馬場様
ありがとうございました。
非常に助かりました
馬場様
初めまして。Gyoと申します。
こちらの過去のコメントにもありましたが、auto.arimaの結果が書籍と大きく異なります。
forecastのバージョンもすでに8.4ですが相変わらずのようです。
trace = Tで確認しましたが、(p,d,q)のdが常に0、(P,D,Q)のDが常に1になります。
これはバグなのでしょうか？それともd, Dはこちらで何かしら設定する必要があるのでしょうか？
stepwise = F としてもこのような結果になるのであれば、auto.arimaを今後も使用して良いものか悩んでしまいます。
以上、お手数ですがご回答いただけましたら幸いです。
よろしくお願いいたします。
Gyo様
コメントありがとうございます。
管理人の馬場です。
どうやらforecastパッケージのauto.arimaのアルゴリズムそのものも多少変わったようです。
書籍と同じ結果を出す場合は、お手数ですがパッケージのバージョンを書籍執筆時の8.2にして実行されるようお願いいたします。
ご回答ありがとうございます。また、返信が遅くなり申し訳ありません。
バージョンを8.2に戻そうと上の方に書かれている通りに試したのですがうまくいきませんでした。
個人的には書籍と同じ結果は出なくても考え方が正しいのであれば構わないと思っていますが、仮に今のauto.arimaにバグがあるのであればそもそも使い続けるのは難しいと感じた次第です。
馬場様
はじめまして，Yudaiと申します。
時系列隼本のVARモデルの章を読んで，自分が持っているデータにもVARモデルを使ってみようかと思ったのですが，データの前処理の段階で足踏みをしています。
書籍を読んでもいまいちピンと来なかったのですが，多変量時系列データにVARモデルを適用しようとしたときは，元データをどのような性質を持つデータに変換するとよいのでしょうか？
Yudai様
コメントありがとうございます。
返信が遅れてすいません。
『元データをどのような性質を持つデータに』
これは難しいところですが、線形かつ定常であれば、
安全かと思います。
（正確に言うならば「VARモデルが満たすべき仮定をすべて持つデータ」となるのでしょうが、
それではあまりに不親切なので……）
とはいえ、VARモデルを適用するだけならば、
単位根があっても（非定常でも）問題ありません。
単位根を持つデータに対してVARを適用して、
その係数のt検定を行っても問題ないということです。
ただし、Grangerの因果性検定を行う場合は、
定常性を満たす必要があります。
というわけで、以下のような前処理になるでしょうか。
1．データをみて、必要なら対数変換(p39参照)
2．Grangerの因果性検定をするならば、Box-Jenkins法と同様に単位根を調べて、
必要ならば差分をとる。
あくまでも「目安」程度ですが。
馬場様
お世話になっております。
以前,質問させていただきましたKです。
auto.arima関数の出力結果に関して不明な点がございます。
auto.arima関数で得られた結果からARIMAモデルを定式化したいと考えたときに,
定数項はどこで確認すればよいかわかりません。
お手数ですがご教授いただけますと幸いです。
馬場様
お世話になっております。
以前、質問させていただきましたKです。
auto.arima関数について質問があります。
auto.arimaの出力結果で定数項はどこでわかるのでしょうか？
お手数ですがご解答のほどよろしくお願い致します。
K様
コメントありがとうございます。
管理人の馬場です。
auto.arimaは、デフォルトでは平均項を入れるか入れないかも含めて、AICで判断します。
想像するに、AIC最小モデルにおいて、平均項が不要になったのではないかと思います。
例えば以下のコードを実行すると「mean」というパラメタが得られるはずです。
library(forecast);
set.seed(1);
ts.sim <- arima.sim(list(order = c(1,0,0), ar = 0.7), n = 200) + 5; mod <- auto.arima(ts.sim, trace = T); mod ちなみに、ts.simの末尾の「＋５」をなくすと、「mean」というパラメタがなくなるはずです。
馬場様
迅速なご返事ありがとうございます。
イメージでは ARIMA with constant
などで出力されると思いましたがmeanで
出力されるんですね。
私が分析しているデータの値が大体60～90だったので、ARMAモデルの平均の式から定数項がないのは何かおかしい気がすると考えましたが、何か私の理解で間違えてる部分があるのでしょうか…。
また、forecastでの予測は最適予測ですか？
もしそうであれば当てはめデータに対して、ホワイトノイズをランダムで発生させながら逐次的に予測する関数などはございますでしょうか？
質問が多くなってしまい申し訳ありません。
お手数ですが何卒よろしくお願い致します。
K様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて失礼しました。
> 私が分析しているデータの値が大体60～90だった
和分過程ですと、一度差分をとるので、データの平均値が大きかったとしても、定数項が入るとは限りません。
> また、forecastでの予測は最適予測ですか？
DGPを完全に推定できているならば、最適予測になるでしょうが、現実的にはこれは達成できていないとみなすところでしょうね。
馬場様
ご返事ありがとうございます。
予測に関して理解がまだ足りないようです。
また調べてみます。
2018/9/27 第3刷で学習しております。
p.136の
prais.winsten(y_ar ~ x_ar, data=d, iter = 1)
を実行すると、
> prais.winsten(y_ar ~ x_ar, data=d, iter = 1)
Error in prais.winsten(y_ar ~ x_ar, data = d, iter = 1) :
could not find function “prais.winsten”
となります。
mod_gls_PW mod_gls_PW
で
Call:
prais_winsten(formula = y_ar ~ x_ar, data = d)
Coefficients:
(Intercept) x_ar
0.4691 -0.0162
AR(1) coefficient rho: 0.814
だけしか出力されません。
praisのバージョンは以下です。
> packageVersion(“prais”)
[1] ‘1.1.0’
何か分かりましたら宜しくお願いします。
安芸様
コメントありがとうございます。
管理人の馬場です。
praisのバージョンが上がったことにより、関数名称や関数の挙動が変わったようです。
以下のコードに変更すると、書籍とほぼ同様の結果が得られるはずです。
mod_gls_PW <- prais_winsten(y_ar ~ x_ar, data=d, max_iter = 1)
summary(mod_gls_PW)
サポートページにも、この件を追記いたしました。
ご連絡いただき、誠にありがとうございます。
有難うございます。
以下のように結果が得られました。
後、馬場さんの本を8日間で読破する事が出来ました。全てのコードも実行しました。齢63になりますが、専門書をこんなに速く読めたのは生まれて初めてです。本当に感謝しております。
> mod_gls_PW summary(mod_gls_PW)
Call:
prais_winsten(formula = y_ar ~ x_ar, data = d, max_iter = 1)
Residuals:
Min 1Q Median 3Q Max
-4.6021 -1.3923 0.0357 1.1703 4.8387
AR(1) coefficient rho after 1 Iterations: 0.8148
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.46450 0.27551 1.686 0.0926 .
x_ar 0.01755 0.04984 0.352 0.7249
—
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.031 on 398 degrees of freedom
Multiple R-squared: 0.0005141, Adjusted R-squared: -0.001997
F-statistic: 0.2047 on 1 and 398 DF, p-value: 0.6512
Durbin-Watson statistic (original): 0.3738
Durbin-Watson statistic (transformed): 1.955
済みません。コピペをミスって頭が抜けてました。
> mod_gls_PW summary(mod_gls_PW)
Call:
prais_winsten(formula = y_ar ~ x_ar, data = d, max_iter = 1)
…
と続きます。
安芸様
コメントありがとうございます。
管理人の馬場です。
そのようにおっしゃっていただけると大変に励みになります。
ありがとうございます。
当方の書籍が少しでもお役に立てたのであれば幸いです。
馬場さま
書籍を購入し、勉強させていただいております。
その中で混乱してしまっているところがございましてよろしければご教示ください。
５部３章のカルマンフィルタに関する部分です。
この章では４種類の誤差がでてきます。
・ 過程誤差
・ 観測誤差
・ 状態の予測誤差
・ 観測値の予測誤差
です。ここで
誤差 ⇒ 真の値と実際の値の差
残差 ⇒ 実際の値と予測された値の差
という理解のもとで本を読み進めていったのですが、３章あたりから
・ 状態の予測誤差 と 過程誤差の違い
・ 観測値の予測誤差 と 観測誤差の違い
がよくわからなくなってしまいました。
例えば、２章１９４ページのローカルレベルモデルの説明図（過程誤差・観測誤差）のあたりは
ふむふむという感じだったのですが
状態の予測誤差・観測値の予測誤差がこの図の中であらわされるのか、
そうだとしたらだとどうなるんだ？？？というようなことでモヤモヤしております。
ご教示いただければと思います。
よろしくお願いいたします。
Kaz様
コメントありがとうございます。
管理人の馬場です。
『誤差 ⇒ 真の値と実際の値の差
残差 ⇒ 実際の値と予測された値の差』
確かに、基本的には上記の理解でよいですが、
状態空間モデルの「過程誤差」に関しては、
少し異なる解釈となります。
「過程誤差」は「状態がどれだけ大きく変化することを許すか」
を決めるパラメータだと考えるとわかりよいかと思います。
「観測誤差」もパラメータとなりますが、
こっちは誤差のイメージ通り「状態と観測値のずれの大きさ」を
表すものとなります。
『章１９４ページのローカルレベルモデルの説明図』でいうと、
四角い状態が、時間が変わって次の状態に移るとき、
どれくらい大きく変化できるかを決めるパラメータが過程誤差となります。
「過程誤差」が小さい、極端な場合ゼロだと、四角い状態は全く上下に動きません。
『状態の予測誤差・観測値の予測誤差がこの図の中であらわされるのか』
という指摘に対しては「あらわされていない」が回答となります。
P194の概念図はあくまでも「状態空間モデルの概念図」です。
「状態の予測誤差・観測値の予測誤差」は、状態空間モデルを使って、
実際に、データから状態を推定するときに、
考える概念なのだと思うと良いかもしれません。
ただ、この辺りはどうしても日本語の限界で、
P208の数式において、σ_wとσ_v(過程誤差の大きさ・観測誤差の大きさ)
そしてP_tとF_t(状態の予測誤差の大きさ・観測値の予測誤差の大きさ)
で把握したほうが間違いがないかと思います。
ローカルレベルモデルの場合、
σ_wとσ_vは時間によって変わらないパラメータだというのが、
数式を見るとよりはっきりしますね。
・・・日本語で頑張ると
●「過程誤差」のみ、誤差のイメージから外れる。
●「過程誤差」「観測誤差」はともにパラメータであり、
P194の概念図にも載っている。
●「状態の予測誤差」「観測値の予測誤差」は、
状態を推定・予測したり、観測値を予測したりするときに
現れるものだ。
というのが、日本語での解釈(分類)になるでしょうか。
馬場さま
ご回答ありがとうございました。
１９４ページのイメージの中では状態の予測誤差・観測値の予測誤差については表現されていない
ということで理解しました。
「過程誤差」についてもご説明からイメージを理解しました。
また、「観測誤差」については
誤差を「真の値と予測値の差（⇒単位が同じものの差）」というとらえ方でみていたのですが
ご説明の内容をよんで、観測誤差は「状態（例えば購買力）と観測値（例えば売上金額）の差」
ということだから、状態と観測値の両者に関連はあってもそもそも同じ現象・事象とは限らないという
認識をあらためてしました。 ⇒ 認識がおかしければご指摘ください。
一方、「状態の予測誤差」「観測値の予測誤差」ですが
＞「状態の予測誤差」「観測値の予測誤差」は、
＞状態を推定・予測したり、観測値を予測したりするときに
＞現れるものだ。
とご回答いただいているのですが、「予測【誤差】」が何と何の差なのかまだピンときておりません。
書籍では
＞「観測値の予測【残差】 ＝ 実際の観測値 － 観測値の予測値（２０６ページ）」
とあり、これは納得しております。
一方、こちらのホームページ内の「予測の評価方法：誤差の指標とナイーブな予測」という記事の中には、
＞予測誤差はシンプルに
＞「実際のデータ」 － 「予測値」
＞から計算ができます。
という記述があります。
これは前述（観測値の予測残差）の【残差】に相当するものではないかと思い、
予測誤差と予測残差の違いで混乱しております
（とくに予測誤差が何と何の差なのか）。
同じ「予測誤差」でも「観測値の予測誤差」と「状態の予測誤差」で「予測誤差」のイメージが
違うようでしたらそれぞれ教えていただければ幸いです。
よろしくお願いいたします。
Kaz様
コメントありがとうございます。
管理人の馬場です。
> 状態と観測値の両者に関連はあってもそもそも同じ現象・事象とは限らないという
> 認識をあらためてしました。 ⇒ 認識がおかしければご指摘ください。
これは難しいところですが、若干違うような気がします。
状態は「潜在変数」です。
すなわち、私たちはそれを観測することができないが、
それが存在するとみなそうじゃないか、
ということで登場するものです。
潜在変数というのは、私たち人間が、
現象を理解・解釈しやすくするために導入するものです。
「センサーで記録されたナイル川流量」が観測値なら、
目に見えない「本当のナイル川流量」を状態とみなすような感じですね。
両者は異なるように思えますが、
まったく違うものかといわれると、そうではないはずです。
目に見えない状態と、観測できる観測値とのずれの大きさを、
観測誤差というパラメータで表現したということです。
「状態は観測できない」という事実、
「状態は私たちが”あること”とみなしたものだ」
という事実を少し意識しておくと、
残差という言葉がなかなか出てこない理由が見えてくるかと思います。
(場合によっては状態が本当に”ある”こともあるでしょう。
制御なんかだとセンサーデータから正しい状況を知るわけだから、
状態が本当に”ある”こともあるわけです。
でも、それが計測できないわけだから、
ここは”あることとみなす”と記しても間違いではないはずです。
統計学の文脈だと、本当に”あることとみなされた”というような
どうあがいても計測できなさそうな状態も出てきます)
> 「予測【誤差】」が何と何の差なのかまだピンときておりません。
予測誤差に関しても、状態を考える場合は、
「状態は観測できない」わけですから、
観測値との差分を得ることはできません。
この辺りが残差との違いです。
(そして日本語の限界といえるかもしれない。
誤差という日本語のイメージからみると、確かに少し変かもしれませんね)
p208の数式に戻ると、
「状態の予測誤差の分散P_t」と
「観測値の予測誤差の分散F_t」
と記載してあるはずです。
”目に見えない状態”を相手にして残差は計算できません。
”理論上の予測誤差の大きさ”を表現するのに、
P_tとF_tの2つが登場しているわけです。
誤解を恐れずに言うと、こいつらは
「私たちが”あること”と仮定した存在」
といえるかと思います。
じゃあそんな「私たちが勝手に”あること”とみなした潜在変数(状態)やら○○誤差」
なんてものをどうやってデータと対応付けるのか。
そこでようやく、最後に「残差」が登場するイメージですね。
観測値の予測値は(状態方程式と観測方程式から)計算できる。
なので、観測値の残差も計算できる。
この残差の「分散の大きさ」は、F_tで表現できるはずだ。
ここでようやく、私たちが(勝手に？)想定した、
○○誤差の分散やらが、
データと対応付けられるわけです。
日本語だとどうしても正確に理解するのが少し難しいかもしれませんが、
イメージがつかめれば幸いです。
馬場様
突然のご連絡ならび、他の質問者様へのご回答への追加のご質問、誠に恐れ入ります。私も書籍を購入しまして、現在勉強中なのですが、こちらのスレッドの質問者様と同様に「観測値の予測誤差」と「観測値の予測残差」の違いが分からず、色々調べる中でこちらにたどり着き、質問させていただいた次第です。
スレッドの内容は読ませていただいたのですが、まだ理解しきれておらず、お手すきの際にご回答いただけると幸いです。
まず、「観測値の予測残差」と「観測値の予測誤差」は別のものとしてして認識しており、
「観測値の予測残差＝実際の観測値-予測した観測値」
である一方で、
「観測値の予測誤差＝補正前の状態の予測値-予測した観測値」
と考えたのですが、これは誤りでしょうか？
もし誤りである場合は、式で書くと何と何の誤差なのでしょうか？
また、「状態の予測誤差＝補正前の状態の予測値 – 1期前のフィルタ化推定量（補正後の状態）」と考えたのですが、こちらも認識に齟齬があるでしょうか？
p.220では、「観測値の予測”残差”」の分散が「観測値の予測”誤差”の分散」となっているので、ここでは誤差と残差の分散は一緒であるのか、混乱してしまいました。
取り留めのない質問で申し訳ございません。
ご回答いただけますと幸いです。どうぞよろしくお願いいたします。
AT様
コメントありがとうございます。
返信が遅れて失礼いたしました。
■1つ目の質問：「観測値の予測残差」と「観測値の予測誤差」
誤差と残差の区別は難しいですね。
残差の理解はご指摘の通り「観測値の予測残差＝実際の観測値-予測した観測値」です。
一方の誤差は「実際の観測値-真の値」くらいの意味です。真の値と言われても、そんなもの観測できませんので、誤差は観測できません。
そのため、誤差の代わりに残差を使います。
誤差は理論上の値、残差は実際に観測された値、くらいに考えるとわかりやすいかなと思います。
■2つ目の質問：状態の予測誤差
状態は観測不可能なので、状態の誤差も観測は不可能ですね。
けれども、「”理論上の”状態の予測誤差の分散」は計算できます。
誤差を観測できないのに、誤差の分散の大きさだけは計算できるというのは不思議な気持ちですが、分散はモデルから導かれるのだと考えるといいかなと思います。
そして「”理論上の”状態の予測誤差の分散」をもとにして「”理論上の”観測値の予測誤差の分散」が計算できます（p206）。
そして「”理論上の”観測値の予測誤差の分散」は、「”観察可能な”観測値の予測残差の分散」と等しいはずです。
これを利用してp220のように最尤法を実行しています。
参考になれば幸いです。
馬場さんの本は全て購入して読んでいます。
非常に分かりやすくて勉強になっています。
状態空間モデルの実装について質問があります。
ローカルレベルモデル等を使うときに、
同じモデルから生成された時系列が複数固ある場合、
これら全てを使ってモデルの調整をしてくれる方法は
あるでしょうか。
そのようなパッケージを探しているのですが、見つかりません。
自作するしかないですか。
SM様
コメントありがとうございます。
管理人の馬場です。
書籍をお読みいただき、ありがとうございます。
専用のパッケージといわれると難しいですが、
例えば以下の記事ではStanを使って
「複数の観測値があるが、状態は1つしかないモデル」
を推定しています。
参考になれば幸いです。
「Stanで推定する多変量時系列モデル」
馬場様
お返事ありがとうございます。
こちらの説明不足でした。
観測値が多変量ではなく、
例えば、２つの時系列、
x1=[x_11,x_12,…,x_1n]
x2=[x_21,x_22,…,x_2m]
がある場合です。
この場合、[x1,x2]のリスト等をモデル調整
に渡したいのです。
x1とx2を結合してしまうと、短い時系列のときは
駄目な気がします。
このような複数の時系列を持っているときでも
EMアルゴリズムは適用できると思うので、状態空間モデル
のパッケージでその様な仕様のものがありそうな気がしていました。
SM様
コメントありがとうございます。
管理人の馬場です。
申し訳ないですが、対象となるデータのイメージがつかめませんでした。
２つの時系列がある、というのは、多変量ということではないのでしょうか。
可能なら、データの具体的な名称など教えてください。
また、「モデル調整に渡す」という言葉もあまり一般的ではないかと思います。
こちらは「モデルを推定する」くらいの意味でしょうか。
質問を質問で返してしまって恐縮ですが、
よろしくお願いいたします。
馬場様
度々すみません。
畑違いのため、用語が正確ではないことをお許し下さい。
例えば、気温の１日毎の時系列データがあり、所々途切れている場合です。
２０日間と３０日間の２つの時系列があり、
それぞれ、同じモデルとパラメータから出力されているとします。
この場合、２０日間と３０日間の両方の時系列を観測変数として、
モデルパラメータを調節したいのです。
つなぎ合わせて５０日間のデータとするのが簡単ですが、
これだと接合部分がおかしいですよね。
短い時系列データが沢山あったときに、全てのデータを使って
ローカルレベルモデルのパラメータ推定をすることは
可能でしょうか。
知識不足でお手数をおかけしていてすみません。
SM様
コメントありがとうございます。
管理人の馬場です。
情報ありがとうございます。
こちらのケースですと、「所々途切れている」部分を欠損値として扱うのが良いかと思います。
例えば、1月1日～20日まで観測値が得られていて、
そのあと、2月1日～2月28日まで観測値が得られていたとします。
間の10日ほどが存在しないわけですが、それを欠損として扱います。
1月01日～1月20日まで：データあり
1月21日～1月31日まで：欠損
2月01日～2月28日まで：データあり
のような形で、1月と2月の気温のデータを用意してやり、
これを使ってローカルレベルモデルなりのパラメータを推定すればよいかと思います。
馬場さま
早速の回答をありがとうございます。
ご説明いただいた内容がなんとなくイメージでき、モヤモヤがかなり晴れた気がします。
ありがとうございました。
著者のかたに直接質問できるのは大変ありがたいです。
今後のご活躍もお祈りいたします。
馬場様
素早いお返事ありがとうございます。
欠損値とするのもあり得ます。
しかし、がばっと欠けているところがあるのです。
Statsmodelsでは１つの時系列しか入れる方法をしらないので、
他のパッケージではどうかと思いました。
隠れマルコフモデルの場合、
hidden_markovというパッケージが複数の時系列データを入力でき、
さらに重みまでつけることができるようです。
https://hidden-markov.readthedocs.io/en/latest/example.html
隠れ変数が連続変数の場合でも同じことができると思うので、
探していたところでした。今のところなさそうなので、
しばらくは諦めようと思います。
大変ありがとうございました。
馬場さん
連投すみません。
先の私の返信は勘違いしています。
欠損値は多くても大丈夫でしょうか。
また、Statsmodelsのローカルレベルモデルでも
欠損値は使えるでしょうか。
SM様
コメントありがとうございます。
管理人の馬場です。
「大丈夫」かといわれると難しいですが、
モデルを推定すること自体は可能です。
あまりにも間が長くあいていたら、
そもそも同じデータ生成過程から得られたとみなしてよいか、
難しいところもあるかもしれませんね。
もちろん、欠損が多ければ多いほど、推定結果の不確実性は増すはずです。
結果をどのように解釈・活用するかはケースバイケースでしょう。
ちなみにPythonのstatsmodelsでも欠損値としてnumpyの「NaN」を指定すれば、
そのまま補完をしてくれるはずです。
馬場様
返信遅くなりましたが、お手数ありがとうございました。
それを試させていただこうと思います。
本（第1版4刷）を購入させて頂き、勉強させて頂いています。
高校途中で数学は挫折したので詳しい式になるとよくわからないですが、
概要だけはわかりやすく読めている、ような気がします。
さて、112Pの石油平均価格を予測に使うところ（白枠四角内のひとつめ）の一番下のコードを実行すると
xregはmatrixかvectorにしろ（xreg should be a numeric matrix or a numeric vector）と言われて
しまいます。
素人っぽい質問で申し訳ないのですが、変換の方法がよくわからずこちらご教授頂ければと思います。
一応as.vector等は試してみました。
よろしくお願い致します
※既出質問でしたら、過去回答場所（日時）等頂ければコメント内から探させて頂きます。
YS様
管理人の馬場です。
回答が遅れまして、失礼いたしました。
コメントありがとうございます。
拙著をお読みいただき、ありがとうございます。
当該エラーですが、当方でも再現できました。
forecastパッケージのバージョンが上がったのが理由のようです。
以下のようにas.matrix関数を適用することで対処できるはずです。
forecast(sarimax_petro_law, xreg=as.matrix(petro_law_mean))
サポートページにも、この件を追記いたしました。
ご連絡いただき、誠にありがとうございます。
書籍(第1版3刷)を購入させていただき、勉強させていただいております。
P.273について質問させてください。
最新のRのバージョンでは、”Nippon”パッケージは使えないのでしょうか。
Rver: 3.6.0
インストール時（メッセージ）：
package ‘Nippon’ is not available (for R version 3.6.0
使えない場合、なにか別の方法があれば教えていただけますか。
何もかも初心者で申し訳ありませんが、ご回答よろしくお願いいたします。
たぬき様
コメントありがとうございます。
管理人の馬場です。
当方でも同様の現象を再現できました。
どうやらNipponパッケージがCRANからなくなっており、もはやインストールができなくなってしまったようです。
Nipponパッケージを使わないで同様の処理を行う方針になります。
以下、コードが続いて恐縮ですが、『内閣府：「国民の祝日」について』で提供されている、祝日のCSVファイルを用いる方法を載せます。
p273の2つ目のコードから以下に差し替えれば、書籍と同じ結果が得られるはずです。
# 祝日の取得
holidays <- read.csv("https://www8.cao.go.jp/chosei/shukujitsu/syukujitsu.csv")
colnames(holidays) <- c("date", "holyday_name")
holidays$date <- as.Date(holidays$date)
head(holidays)
# 祝日か否かの判断
# p273の2つ目のコードに相当
head(dates %in% holidays$date)
# 3月21日は日曜日で、3月22日が振り替え休日
dates[dates %in% holidays$date]
# 曜日
weekdays(dates[dates %in% holidays$date], T)
# 日曜日は取り除く
holiday_date <-
dates[(dates %in% holidays$date) & (weekdays(dates, T) != "日")]
holiday_date
# 祝日フラグ
holiday_flg <- as.numeric(dates %in% holiday_date)
holiday_flg
後ほど、別の記事などにおいて、祝日判定の方法を解説いたします。
まずはこちらのコードを参照いただければと思います。
馬場様
隼時系列本（第1版第4刷）を拝読し、状態空間モデルについて勉強しています。
他の時系列系の書籍と比べても大変理解しやすく、重宝しています。ありがとうございます。
宜しければ1つ質問させて下さい。
KFASを使って何らかのモデルを構築し、まず手元にあるデータで学習をしたとします。
その後新しいデータを入手したので、上記学習済みモデルに追加学習させたいとします。
このような場合、どのようにコーディングしたらよろしいでしょうか。
ネットで調べても追加学習についての解説は見つけることができませんでした。
御本の範囲を超えた質問で恐縮ですが、可能な範囲でご教示いただけますと幸いです。
宜しくお願い致します。
Yoshiki様
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただき、ありがとうございます。
追加学習というのは、深層学習など、
学習に時間がかかるものに対して適用されることが多いかと思います。
KFASでのパラメータ推定は、そこまでコストがかかるものでもありませんので、
あたらしいデータを用いて、パラメータの推定を最初からやり直してもよいかと思います。
馬場様
やはりKFASには追加学習するという機能はないのですね。
状態空間モデルの原理的には可能だという理解なのですが、
確かに最初からモデリングし直すという方法もありますね。
ご回答、誠にありがとうございました。
本書を使って時系列モデルを初めて学習しました。、
大変に読みやすく、応用・実装例も豊富で参考になっています。
１点だけ、気になったのですが、p289の「2-6 補足：確率の基本公式」の周辺確率・同時確率・条件付確率の説明の部分で、P(θ=女性|y=赤鞄発見)= … となっていますが、説明の流れ的には、P(θ=赤鞄所持|y=女性)= … ではないでしょうか？
第４刷を参照しております。
toshi様
コメントありがとうございます。
管理人の馬場です。
ご指摘ありがとうございます。
確かに、この一文は、説明の流れとして不整合がありそうです。
P(y=赤鞄所持|θ=女性)= …でしょうか。
出版社と調整のうえ、修正を検討いたします。
誤植がありまして、誠に申し訳ございません。
ご指摘いただけたこと、感謝いたします。
当サイトとご著書に大変お世話になっている者です。「共和分」についてわかり易い内容の本で勉強したいと思っています。「時系列分析と状態空間モデルの基礎」が良いのではと考えましたが目次には無く、書名と共和分で検索しましたが、見当たらないようです。この本での「共和分」についてのご解説の有無について教えて頂けるととても助かります。よろしくお願いいたします。
うどん様
コメントありがとうございます。
管理人の馬場です。
共和分に関しては、あまり量はありませんが、
「第3部第1章：見せかけの回帰とその対策」で、5ページほど簡単な紹介を載せています。
https://logics-of-blue.com/wp-content/uploads/2018/01/tsa-ssm-book-contents.pdf
馬場 様
”共和分”についても詳細な目次には表示されていたのですね。教えて頂きありがとうございました。ご著書を利用させて頂きます。いつもありがとうございます。
隼時系列本でベイズ推論とStanについて学んでいます。とてもわかりやすく、Rで実践もできて勉強になっています。ありがとうございます。
恐れ入りますが1つ質問をさせて下さい。
第6部4章で、最後に観測誤差を除いた場合の個体数の変動のみを図示されていますが、観測誤差を含めた場合を見ないのはなぜでしょうか。
人による観測は誤差が大きいので、真の個体数の変動を推定するために観測誤差を除いた場合の期待値を見ている、という理解で正しかったでしょうか。
※実際にやってみたところ、観測誤差を含めるとほぼ観測値と一致してしまいました。これを解釈するため第5章213ページの下部の図を参照すると、観測誤差ではなく過程誤差が大きいことになりそうで、少し混乱しています。
どうぞ宜しくお願い致します。
Yoshiki様
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただきありがとうございます。
1．観測誤差を除いた意図について
観測誤差を除いて、全体の傾向(トレンド)を確認するのが大きな目的でした。
2．以下コメントへの返信
> ※実際にやってみたところ、観測誤差を含めるとほぼ観測値と一致してしまいました。
> これを解釈するため第5章213ページの下部の図を参照すると、
> 観測誤差ではなく過程誤差が大きいことになりそうで、少し混乱しています。
p213の図は、状態の推定値を太線で示しています。
状態の推定値は「観測誤差を除いた値」とみなしても良いでしょう。
p213の下図で、太線と細い線が一致しているのは「観測誤差を除いたはずの状態推定値が観測値と一致している」ということです。
そのためp213の下図は過程誤差が大きくて、観測誤差は小さいことになります。
一方の第6部第4章のモデルにおいて、
観測誤差を除いて、全体の傾向(トレンド)を確認した図(書籍p328に載っている図)では、
観測値と状態推定値が一致していません。
なので、観測誤差が大きいことになります。
ここに矛盾はないものと考えますが、いかがでしょう。
馬場様、返信ありがとうございました。
ご解説を読んで、私が誤って理解している部分に気づきました。改めて書籍も読んで理解いたしました。
ご解説、ありがとうございました！
隼時系列本を読ませていただいております。
今まで未知の領域であったカルマンフィルターが少しずつ頭に入ってきており
本当に感謝でいっぱいです。
私がやりたいことは、毎日の気温の予想です。
気象庁の予報データと実際の観測値にはずれがあるため
その補正をしたいと考えています。
気象庁の資料をみると、「線形ガウス状態空間モデル」、「ベイズ推定」
が今後のキーワードになると思うのですが、実装面も含めて
次に読むべき本を指南いただけると助かります。
どうぞよろしくお願いいたします。
馬場様
時系列分析と状態空間のモデルの基礎ですが、大変勉強になっております。
次に読むべき書籍の指南をお願いしたいと考えています。
具体的にはカルマンフィルターを用いて、気温の予報（気象ガイダンス）を
したいと考えております。
気象庁の資料をみると、『線形ガウス状態空間モデル』、『ベイズ推定』といった
内容がキーワードになると思います。
どうぞよろしくお願いします。
ありた様
コメントありがとうございます。
また、拙著をお読みいただき、ありがとうございます。
2件ご質問がありましたが、おそらく同じ内容と思われるので、
こちらだけ返信いたします。
気象の予測に関しては、当方は行ったことはありません。
ありた様の課題と対応しているかはわかりませんが、
例えば「岩波データサイエンスVol.6」では、
時系列分析と気象におけるデータサイエンスの話題が載っています。
気象というテーマをぬいて、純粋な状態空間モデルの教科書でいえば
『基礎からわかる時系列分析』
https://github.com/hagijyun/tsbook
が難易度的に、拙著の次に読む本としてお勧めできます。
参考になれば幸いです。
馬場様
隼本大変わかりやすく、楽しみながら読ませていただいております。
質問なのですが、第二部2-7対数変換とその解釈の部分で、
時系列データ = 周期変動 + トレンド + ホワイトノイズ
と考え、これを対数変換すると
log時系列データ = log 周期変動 + logトレンド + log ホワイトノイズ
と書かれています。
これは、数学的にはおかしい気がしますが、意味としてはどちらかといえば
log時系列データ = 新周期変動 + 新トレンド + 新ホワイトノイズ
という感じなのでしょうか？
もしお手すきでしたらご教授いただけると幸いです。
TN様
コメントありがとうございます。
管理人の馬場です。
返信が遅れてしまい、大変失礼いたしました。
こちらですが、ご指摘の通り、記載の仕方がぶれているように思います。
申し訳ございません。
こちら、記載の仕方を検討いたします。
ご連絡いただきまして、誠にありがとうございました。
何度も読み返して確認しています。
3-3 ARCH-GARCHモデルで以下発生しました。以前は稼働していました。
2021/4/9 Error in library(rugarch) :
‘rugarch’ という名前のパッケージはありません。
が出てライブラリが起動しなくなりました。
再度インストールしましたが変化ありません。
対応策はありますか？
藤田様
コメントありがとうございます。
管理人の馬場です。
返信が遅れて失礼いたしました。
当方の環境ではライブラリの読み込みができました。
古いバージョンのRを使っていると、読み込みに失敗するかもしれません。
新しいバージョンのR(例えば4.0.5など)の利用を検討してみて下さい。
はじめまして、
現在カルマンフィルタのシステムノイズがガウス分布ではなくコーシー分布であるものの実装を試みておりますが、拡張カルマンフィルタで解くのかなと思いつつ、具体的な手順を記載している書籍等がなく困っております。
こちらの書籍にそういった内容はございますでしょうか。
教えていただければ幸いです。
よろしくお願いいたします。
shimatani様
コメントありがとうございます。
返信が遅れて失礼いたしました。
拙著には拡張カルマンフィルタの解説は載っていません。
参考になれば幸いです。
初めまして,
本質的ではない質問で恐縮ですが, 気になったので投稿させていただきました.
p.39の最後の式,
時系列データ = 周期的変動 + トレンド + ホワイトノイズ
の部分ですが,
次ページで「データを対数変換」とあるので, 両辺に対数を取ったと考えたのですが,
log時系列データ = log周期的変動 + logトレンド + logホワイトノイズ
となっていました.
これはどういう式変換をしているのでしょうか？
よろしくお願いいたします.
momomo様
コメントありがとうございます。
管理人の馬場です。
拙著をお読みいただきありがとうございます。
ご指摘の箇所、数式の上ではこれで正しいのですが、
すこし分かり難い解説になってしまったかもしれません。
ご指摘の通り、通常は、データに対して対数変換をする場合
「周期的変動」や「トレンド」「ホワイトノイズ」を個別に取得はできません。
ここでは、「仮にこれらが個別に取得できた場合、各々の対数変換をすると、掛け算とみなせる」くらいの解釈になります。
あくまでも、モデルの解釈上の話であることに注意してください。
逆向きの解説をした方がわかりやすかったかもしれません。
すなわち、もともとのデータが
時系列データ＝周期的変動×トレンド×ホワイトノイズ
という掛け算のモデルで表現できたとします。
ここで対数変換をすると以下です。
log(時系列データ)＝log(周期的変動×トレンド×ホワイトノイズ)
log時系列データ = log周期的変動 + logトレンド + logホワイトノイズ
要するに、対数変換をしたデータlog(時系列データ)に、加法的（足し算）な構造のモデルを適用した場合、
それは乗法的（掛け算）な構造のモデルとみなせるということです。
このとき、加法的な構造における個別の項には対数がかかります。
馬場 様
早速ご返信いただきありがとうございます.
なるほど, そういう意味だったのですね.
丁寧に解析いただきありがとうございます.
理解ができました.
馬場 様
お世話になっております.
p.136のPrais-Winsten法の, Rによるpraisパッケージを用いた実行の件です.
2019年2月26日に注意事項で, prais_1.1.0からprais.winsten関数がprais_winsten関数に変更になったということで,
> mod_gls_PW = prais_winsten(
y_ar ~ x_ar,
data = d,
max_iter = 1
)
> summary(mod_gls_PW)
を実行したのですが,
prais_winsten(y_ar ~ x_ar, data = d, iter = 1) でエラー:
引数 “index” がありませんし、省略時既定値もありません
となってしまいます.
> ?prais_winsten
をしてhelpを見て,
> d = data.frame(
y_ar = y_ar, x_ar = x_ar,
time = 1:n_sample # 追加
)
> mod_gls_PW = prais_winsten(
y_ar ~ x_ar,
data = d,
index = “time”, # 追加
max_iter = 1
)
> summary(mod_gls_PW)
としたのですが, この出力では, p.137の
[[2]]
Rho Rho.t.statistic Iterations
0.8078253 27.14816 1
が得られませんでした.
何か対応策はありますでしょうか？
よろしくお願いいたします.
momomo様
コメントありがとうございます。
管理人の馬場です。
コードが動かない件、情報ありがとうございます。
また、実行方法についてもご記載いただき、ありがとうございます。
頂いたコード、可能ならブログ本文にも記載して良いでしょうか？
> 何か対応策はありますでしょうか？
ざっと調べたところ、まったく同じ結果を出す方法は無いようです。
おそらく気になるのはRhoの値かと思います。これに関しては「mod_gls_PW$rho」で取得できるようです。
馬場 様
ご返信いただきありがとうございます.
下記対応策, ありがとうございます.
Rhoはそのようにして出せたのですね.
コードの件, 載せていただいても問題ありません.
少しでもお力添えできれば幸いです.
馬場様
初心者です。お世話になります。
216ページの（5-24）の、ｒt-1= の式がどのように導出されたのかがよく判りません。
214ページ～215ページの「日本語で読む平滑化」を読んでも、この部分の説明は無いように思います。
同様に（5-25）の、ｒT-1= の式展開で、左辺=中央の式=右辺 の、中央の式=右辺 になる理由もよく判りません。
中央の式には（1-KT）ｒT の項が加わっており、右辺には（1-KT）ｒT の項が無くなっています。
にもかかわらず両者が等しい、ということは（1-KT）ｒT = 0 ということでしょうか？
多分、式の意味を私がよくわかっていないせいだと思うので、（5-24）（5-25）の式の解説をよろしくお願いいたします。
ihara_e様
コメントありがとうございます。
管理人の馬場です。
> 「日本語で読む平滑化」を読んでも、この部分の説明は無いように思います
この点はすいません。日本語での説明にはやや限界がありまして、ここの説明は端折りました。
スムージングの導出はやや難解でして、端折る教科書も多いと思います。
平滑化漸化式については「野村(2016)カルマンフィルタ」が、平滑化の導出については「森平(2019)カルマンフィルター入門」などが参考になるかと思います。
> 同様に（5-25）の、ｒT-1= の式展開で、左辺=中央の式=右辺 の、中央の式=右辺 になる理由もよく判りません。
これは簡単で、式(5-24)の直下にかかれているように、漸化式のr_Tの初期値を0と定めているからです。r_Tが0なので、（1-K_T）r_Tも0になります。
馬場様
（5-24）につきまして、参考文献をご連絡いただき、ありがとうございました。
入手して、該当部分を勉強してみます。
また（5-25）につきまして、うっかりしていてすみませんでした。確かにr_T=0と書いておられますね。
とても助かりました。
馬場 様
お世話になっております.
p.216とp.217の状態平滑化漸化式のところで質問があります.
(5.24)式の2行目では,
\hat{\mu}_t = \mu_{t|t} + P_{t|t} r_t
となっていたので,
(5.27)式の1行目は, (5.24)式でt=T-2を代入し,
\hat{\mu}_{T-2} = \mu_{T-2|T-2} + P_{T-2|T-2} r_{T-2}
になると考えたのですが,
本書では,
\hat{\mu}_{T-2} = \mu_{T-2|T-2} + P_{T-1|T-2} r_{T-2}
と, Pの添字が, 「T-2」ではなく「T-1」となっていました.
どうして「T-1」なのでしょうか？
よろしくお願いいたします.
馬場 様
お世話になっております.
p.275の「KFASによる基本構造時系列モデル」の実装における推定結果の確認で,
周期成分の推定結果を確認する際に”sea_dummy1″を描写しているかと思います.
sea_dummy1のみを表示しているのは, 他の5つの周期成分のパラメタ推定値は
sea_dummy1の値が対応する時点ずつズレているだけなので, そうしていると思うのですが,
sea_dummy1の第1成分の値は, 何曜日に対応するものなのでしょうか?
p.201で, 「周期的変動のモデル化」について議論されているかと思いますが,
その箇所と対応付けて考えてみると, sea_dummy1の第1成分の値は火曜日
に対応しているものなのかなと考えました.
また, そう考えていくと, sea_dummy1の第1から6成分を足し合わせた値は
月曜日に対応する値なのでしょうか?
複数質問してしまい恐縮ですが, よろしくお願いいたします.