今回のモジュールでは目的変数に与える色々な外部環境が変化したとき、これまでの予測技術を応用することで影響を評価することを試みる。

本当はボックスジェンズ流のARIMA分析を拡張して、複数のデータの依存関係をみながら予測するのが本筋なのであるが、複数のデータになった途端に学ぶべきことが倍増以上になるのが「多変量時系列の予測」である — 単に差をとってトレンドを消せばよいというわけではなくなる。

この授業ではパッケージ“forecast”で可能な範囲にとどめておくことを留意されたい。

library(forecast)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
##  以下のオブジェクトは 'package:base' からマスクされています: 
## 
##      as.Date, as.Date.numeric 
## 
## Loading required package: timeDate
## This is forecast 6.2

以下、第1節では自動車販売台数がリーマン危機や東日本大震災というショックからどれほどの影響を受けたかを評価する。次に第2節ではNHKの朝ドラ『マッサン』が放映されたことで、ドラマの舞台となった余市町を訪れる観光客数にどれほどのプラス効果があったかを分析した。

どちらも文字通りの将来予測ではないが、環境変化の影響を評価する際にも予測の技術が活用できる点を見てほしい。無条件の予測が役に立つこともあれば、説明変数を含めたうえで複数ケースのシミュレーションを行うこともある。どのように使うかが本モジュールの要点である。

[１]自動車販売：金融ショック・災害ショックの影響評価

自動車販売データはモジュール５の事後課題で分析したことがある。再度、このデータを検討してみよう。

data0 <- read.csv(file="Car-Japan-1993Q1-2015Q2.csv")
head(data0)
##   year quarter    sale crisis
## 1 1993       1 3101281      0
## 2 1993       2 2776345      0
## 3 1993       3 2738075      0
## 4 1993       4 2611844      0
## 5 1994       1 2724604      0
## 6 1994       2 2515763      0

今度はデータファイルの中に“crisis”という系列が含まれている。値は2008年第3四半期から1年間、2011年第1四半期から1年間だけが1であり、他は0である。大きな外部ショックが発生して、その影響の大きさをノイズ（＝不規則成分）と解釈するには無理があるとき、あらかじめ（ザックリと）「平常期」と「危機」の区分を設けたわけである。

【確認】

他に時期を区分する方法はなかっただろうか？

まず読み込んだデータに日付をつけて時系列データにする。

sale <- ts(data0$sale,freq=4,start=c(1993,1))
crisis <- ts(data0$crisis,freq=4,start=c(1993,1))

これが今回の「全データ」になる。

リーマンショックの影響

まず一つめの巨大ショックである「リーマン危機」の影響をデータから引き出すことを考えよう。

危機前の予測と実績との違い

リーマン危機の影響評価だが、最も分かりやすい方法は以下の様なやり方であろう。

  1. 危機発生前のデータだけを用いて、危機発生後の（たとえば）2年間（＝8四半期）を予測する。

  2. 上で得た予測値と実績との差を求める。この差が、2年間に現れたリーマン危機の影響と考える。

要するに「リーマン危機がなかりせば販売台数はこうなっていたはずだ」という経路と実際の経路を比較する方法である。このような方法は一般に経済政策や経営戦略変更の効果を定量的に測定するための定石になっている。ある前提（たとえば実際の政策を維持する）を設けた結果を標準ケースと呼び、他の前提（たとえば政策Aに変更する）を置いた結果をたとえばケースAと呼ぶ。そして「ケースA」と「標準ケース」の違いが政策Aの効果である。そんな議論をするわけである。

【補足】

上の方法においては、現実の経路が標準ケース、危機がなかった場合の予測経路がケースAに該当する。

こんな発想で計算してみよう。この場合は説明変数を用いたARIMAモデルを用いる必要はない。ただ、データを切り分けておく必要がある — データ・ハンドリングはデータ解析の勘所である。

train <- window(sale,end=c(2008,2))
test <- window(sale,start=c(2008,3),end=c(2010,2))

学習データはリーマン危機直前である2008年第2四半期まで。テストデータは危機発生後2年間（8四半期）としている。

学習データを用いてテスト期間である危機後2年間を予測する。

kekka0 <- auto.arima(train)
yosoku0 <- forecast(kekka0,h=8)
bnd <- range(test,yosoku0$mean) # 上限、下限を求めておく
plot(yosoku0,ylim=bnd)
lines(test,col="red")

予測結果をそのまま図に描くと危機発生後のデータが縦軸の目盛に収まらなくなってしまう。そこでテスト期間中のデータ、予測値の両方を含めた値の範囲を求め、それを縦軸の目盛範囲に指定している。

図をみるとリーマン危機の影響は明らかであり、発生後の2年間において、実際の販売台数は95%予測区間外にまで急減している。文字通り、これほどの販売急減は「想定外」であったわけであり、「危機(Crisis)」という形容が当てはまっていた。

予測値と実績値との差がリーマン危機の影響だとすれば、計算は以下のように簡単にすむ。最後に2年間トータルでどの程度の影響があったかを求めている。

eikyou <- test-yosoku0$mean
eikyou
##               Qtr1          Qtr2          Qtr3          Qtr4
## 2008                                  548.5168  -454417.1053
## 2009 -1571019.9247 -1204103.3274  -812571.8106  -617833.4327
## 2010  -676525.2521  -648975.6547
sum(eikyou)
## [1] -5984898

リーマン危機発生後2年間で事前の予測より概ね598万台の販売減少となった。

ある特定の環境変化・制度的変化の影響を評価するとき、上の様なやり方は分かりやすい。

金融危機の影響を測定してから予測する

前節の計算から分かったことは、『通常の予測モデルでは危機発生後の推移を予測できなかった』という点である。

ということは、リーマン危機によるショックを単なるノイズ（＝不規則成分）として処理するには無理があるということだ。危機発生後の一定期間は他とは区別される必要がある。それが今回のデータファイルにある“crisis”というデータである。“crisis”は0と1の二値系列で、このようなデータを回帰分析ではダミー変数と呼ばれることが多い。ダミー変数は0と1の二つの値をもつが、それぞれが意味することは「危機ではない平常期」、「危機である」という質的な情報である。

目的変数に影響を与える説明変数を含めてボックス・ジェンキンズ流の予測モデルを推定したい時は、推定コマンドの中でxreg指定する。

本節では次の設問をとりあげよう。

【設問】

危機発生後1年が経過した時点（＝2009年第2四半期）にあるとする。当該時点以降1年間を予測したい。

前に述べたようにリーマン危機後の販売台数は想定外であって、平常期とは明らかに異なった経路をとっている。故に、危機発生後のデータを平常期のデータと混ぜて、一つのARIMAモデルを当てはめるのは不適切である。

以下のようにデータを切り分けておく。

  1. 2009年第2四半期までを学習データとする。

  2. 2009年第3四半期から1年間をテストデータとする。

データを学習用とテスト用に切り分ける時、“sale”だけではなく、説明変数にする“crisis”も同じようにしておく。

train.sale <- window(sale,end=c(2009,2))
test.sale <- window(sale,start=c(2009,3),end=c(2010,2))
train.crisis <- window(crisis,end=c(2009,2))
test.crisis <- window(crisis,start=c(2009,3),end=c(2010,2))

危機発生後1年間を含んだ学習データに基づき、説明変数を含んだARIMAモデルを推定する。

コマンドはauto.arimaでよい。

kekka1 <- auto.arima(train.sale,xreg=train.crisis)
yosoku1 <- forecast(kekka1,h=4,xreg=test.crisis)

学習期間から予測モデルを推定するので、販売台数をその期間に、また「平常期」と「危機」を区分するデータ“crisis”も学習期間に限定しておく。

次に、将来1年（4四半期）を予測するが、予測期間（＝テスト期間）に説明変数がとる値をxregに指定する必要がある。上で指定した“crisis”の値はデータファイル中の値、つまりテスト期間は「平常期」にあたる、そう前提している点に注意せよ。

予測ラインを描いたグラフ画面を残したまま以下のコマンドを入力して学習期間における予測値、テスト期間の事前予測値を書き加えよう。

plot(yosoku1)
names(yosoku1)　# 計算結果一覧を表示
##  [1] "method"    "model"     "level"     "mean"      "lower"    
##  [6] "upper"     "x"         "xname"     "fitted"    "residuals"
lines(yosoku1$fitted,lty=2)
points(test.sale,pch=20,col="red")

コマンド“forecast”の結果“yosoku1”にどんな結果が含まれているかは、names(yosoku1)を入力して結果一覧を表示させると分かる。データ期間中の予測値の名前が“fitted”であると見当がつくので、これを点線で書き加えた。また、テスト期間中の値は元のデータからtest.saleという名前で切り分けた。これを赤い丸で加えている。

図をみると、危機発生後の販売急減をフォローできている。これは説明変数として“xreg”を含めたためである。反面、テスト期間の予測は大きく外れており、過大な予測になっている。

過大予測になってしまった原因としては、テスト期間中のcrisisの値が0、つまり将来1年が「平常期間」であると前提している点が考えられる。データファイルでは、テスト期間におけるcrisisの値を0としているのだが、もしこの値が1であったなら、つまり危機発生後2年目もなお販売台数は「平常期」とは異なると想定すると、販売台数はどう予測するべきか。こんなシミュレーションも必要であろう。

今回の予測モデルではcrisisを説明変数に含めている。故に、crisisの値が1となった時の販売台数を計算できる。それには以下のように、データファイルとは違う値を“crisis”に指定すればよい。

new.test.crisis <- rep(1,4)　# １を4個ならべる
new.test.crisis <- ts(new.test.crisis,freq=4,start=c(2009,3)) # 4個の1に日付をつける
yosoku1.new <- forecast(kekka1,h=4,xreg=new.test.crisis)
plot(yosoku1.new)
points(test.sale,pch=20,col="red")

今度は金融危機の負の効果が予測に織り込まれたので、過大な予測を出すという誤りを避けることがでた — すべて事前のデータだけを使っており、文字通りの将来予測である点に注意。

以上の作業をまとめておこう。

二つの前提で予測計算をした。危機発生後2年目も危機が続くという前提。もう一つは平常に戻るという前提。この二つの仮定それぞれの下で販売台数を予測した。時間が過ぎれば、事後的に実際の経路が明らかになる。二つの予測を実際の値と比べておこう。

sum(test.sale)
## [1] 9456386
sum(yosoku1$mean)
## [1] 11631116
sum(yosoku1.new$mean)
## [1] 9219361

テスト期間を通して、実績は合計で946万台。最初の予測（＝テスト期間は平常期）が1163万台、二番目の予測（＝テスト期間は危機）では922万台。こんな数字になる。

【補足】

予測とは将来に関する作業である。「こうなると予想できる」というシミュレーションにすぎない — というより、将来に向かってはシミュレーションのみが可能である。その計算には様々な前提がつきまとう。中には重要な前提もある。どのような前提をおくかが予測において非常に重要であるのなら、前提を変更して、何通りかのシミュレーションができるような計算方法をとっておくことが大事である。

繰り返しになるが、二つの予測の差を調べておこう。

tigai<-yosoku1$mean-yosoku1.new$mean
tigai
##          Qtr1     Qtr2     Qtr3     Qtr4
## 2009                   602938.8 602938.8
## 2010 602938.8 602938.8
sum(tigai)
## [1] 2411755

もし2009年第3四半期以降1年間においても金融危機の影響が販売の足を引っ張るなら240万台程度のマイナスを覚悟しなければならない。こんな数字になる。もちろんこの違いはリーマン危機という金融不安の影響からもたらされるものである。

【要確認】

説明変数を含めない単純なボックスジェンキンズ法は環境変化に対して相当頑健である点に留意されたい。試みに、説明変数を含めず、危機から1年が経過した時点以降、将来4四半期を予測するとすれば、どんな予測になっていただろうか。上で説明変数を含めて行った計算と予測精度はどの程度落ちる（？）だろうか？

東日本大震災の影響

次の設問に回答せよ。

【設問】

東日本大震災発生後1年が経過したあとの販売回復について見通しをたてよ。

[２]余市観光に対する朝ドラ『マッサン』の放映効果

月次観光客数のデータファイル

2014年10月から2015年3月にかけて余市を舞台とする朝ドラ『マッサン』が放映され、それを契機に余市町、特にニッカ余市蒸溜所を訪れれる観光客が急増した — 放映終了後も観光客数は高い水準を続けている。

本節では余市観光に与えたドラマの放映効果を評価する。そこでまずは人数ベースでどの程度の効果があったかを検討しよう。

余市町の観光客数データ（月次）から編集されたデータファイルが“Yoiti-Nikka-Others-Monthly.csv”である。そのデータを読みこむことから始めよう。

data0 <- read.csv(file="Yoiti-Nikka-Others-Monthly.csv",header=T)
dim(data0)
## [1] 60  6
head(data0)
##   year month  total nikka others massan
## 1 2010     4  45713 13664  32049      0
## 2 2010     5  85468 25365  60103      0
## 3 2010     6  72567 25989  46578      0
## 4 2010     7 175331 37436 137895      0
## 5 2010     8 137363 34273 103090      0
## 6 2010     9 124779 30722  94057      0
tail(data0)
##    year month  total nikka others massan
## 55 2014    10 147168 74306  72862      1
## 56 2014    11 100928 55082  45846      1
## 57 2014    12  45129 23069  22060      1
## 58 2015     1  46692 24053  22639      1
## 59 2015     2  66917 41541  25376      1
## 60 2015     3  88698 57506  31192      1

データファイルの内容から分かるように、データの開始月は2010年4月、最終月は2015年3月である。また、観光客数は合計、内訳としてニッカ蒸溜所、ニッカ以外のその他合計を訪れた観光客数が含まれている。系列“massan”は放映時期を示すダミー変数（質的情報）である。放映されることが多数の人に意識されるようになった2014年4月から実際に放映された2015年3月までの期間は値1をとり、その他の時期は値0になっている。

これらの個別データに日付をつけて時系列データにするのが次の入力である。

yoiti<- ts(data0[,3:6],freq=12,start=c(2010,4))
total<-yoiti[,1]
nikka<-yoiti[,2]
others<-yoiti[,3]
massan<-yoiti[,4]
par(mfrow=c(2,1))
plot(nikka)
plot(others)

ニッカ蒸溜所とニッカ以外の観光地それぞれを訪れた観光客数の推移をみれば、両者の動きの違いは明らかである。そこで、ニッカ蒸溜所とニッカ以外の観光地を分けて検討を進めることにする。

ここで、観光客合計の年度集計値を求めておこう。

朝ドラ「マッサン」の放映をきっかけにして、余市町を訪れる観光客は2014年度に125万人となり前年度比で45.5%増加した（前年度は5.6%増）。この内、ニッカ余市蒸溜所を訪れた観光客数は55万人で、増加率は93.9%の増加とほぼ倍増した（前年度は7.2%増）。

ドラマ放映効果の計算

【設問】

朝ドラ『マッサン』放映の効果を評価したい。そのための「分析作戦」を提案せよ。利用できるデータは。2014年度末までの余市町観光客数データ（月次）である。

（以下、授業当日に紹介）

効果測定の考え方

実際にはドラマ放映後の2014年度に余市町を訪れる観光数の急増が確認されている。とはいえ、このデータにはドラマ放映とは無関係の観光客も含まれているはずである。単に2014年度の前年度比を『ドラマが放映されたのでこうなった』と言うのみでは正確さに欠けるだろう。

ドラマ放映効果は

  1. ドラマが放映されなかった場合の観光客数
  2. ドラマが放映されている実績を踏まえた将来予測値を接続した観光客数

上の二つの数字の差として定義されるはずである。2014年度、15年度2年間について余市町観光k客数の月次データに基づいて計算することにしよう。

ニッカ観光客の事前予測と放映効果

ニッカ観光客数について、時系列解析モデルを当てはめる（bj0,bj1を推定）。

bj0.nikka <- auto.arima(window(nikka,end=c(2014,3)))
bj0.nikka
## Series: window(nikka, end = c(2014, 3)) 
## ARIMA(0,0,0)(0,1,1)[12] with drift         
## 
## Coefficients:
##          sma1    drift
##       -0.4313  78.4149
## s.e.   0.2886  20.8487
## 
## sigma^2 estimated as 3835284:  log likelihood=-325.19
## AIC=656.37   AICc=657.12   BIC=661.12
bj1.nikka <- auto.arima(nikka,xreg=massan)
bj1.nikka
## Series: nikka 
## ARIMA(0,0,2)(0,1,0)[12]                    
## 
## Coefficients:
##          ma1     ma2     massan
##       0.7351  0.4135  19780.064
## s.e.  0.1346  0.1684   3345.277
## 
## sigma^2 estimated as 30742094:  log likelihood=-482.24
## AIC=972.48   AICc=973.41   BIC=979.97

上の“bj0.nikka”は、2014年3月までのデータに基づく予測モデルであり、マッサンの放映効果が現れる以前のデータを反映するものである（と解釈しよう）。

マッサン放映がなかりせば2014年度、15年度の観光客数はどのように推移したであろうか。これはマッサンの放映効果をまったく含んでいないモデル“bj0.nikka”を利用して2014年度、15年度を予測すればよい。マッサン効果を含まない予測値を2014年度、15年度の標準ケースと呼ぶ。

fr0.nikka<-forecast(bj0.nikka,h=24)
plot(fr0.nikka)

fr1.nikka <- forecast(bj1.nikka,h=12,xreg=data.frame(massan=rep(1,12))) # 2015年度もドラマ放映効果はつづく
plot(fr1.nikka)

次に、モデル“bj1.nikka”は2014年度までのデータに基づく2015年度の予測モデルである。この中には、（もちろん）2014年10月から始まった『マッサン』の放映効果が含まれている。つまり、上の予測値“fr1.nikka”はマッサン放映を踏まえた2015年度の将来予測値である。これを放映後ケースと呼ぶ。

【要検討ポイント】

2015年度の予測をたてるには、上のような方法が唯一・最適の方法だろうか？他に計算の進め方、前提の置き方はないだろうか？

ニッカ以外の観光地

次に、ニッカ以外の施設の観光客数についても同様の計算をしておく。bj0.others, bj1.othersを推定。

bj0.others <- auto.arima(window(others,end=c(2014,3)),D=1)
bj1.others <- auto.arima(others,D=1,xreg=massan)
fr0.others <- forecast(bj0.others,h=24)
plot(fr0.others)

fr1.others <- forecast(bj1.others,h=12,xreg=data.frame(massan=rep(1,12)))
plot(fr1.others)

以上の計算から、「ニッカ」、「その他」の両方について以下の系列を得る。

  1. 標準ケース：マッサン放映がなかった場合の2014年4月‐2016年3月予測値

  2. 予測値：マッサン効果が2015年も継続する場合の2015年4月‐2016年3月予測値。この予測値を2014年度中の実績と接続すれば放映を踏まえた2年間の観光客数が得られる。

ドラマ放映効果の計算（月次）

2014年4月から2015年3月までの期間について観光客合計を算出する。

  1. 標準ケースの合計（total）：2014年4月から2016年3月までは“bj0”による「標準ケース」を使用。

  2. 放映後の合計：　2014年4月～2015年3月までは実績、2015年4月～2016年3月まで“bj1”による「放映後ケース」を使用。

放映効果は
放映効果＝放映後ケース－標準ケーｽ

とする。

まず、標準ケースの観光客合計を2014年4月から2016年3月までの期間について求める。

hyoujun.nikka <- fr0.nikka$mean
hyoujun.others <- fr0.others$mean
hyoujun.total <- hyoujun.nikka + hyoujun.others
hyoujun.total
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2014                                38327.96  71299.28  81841.06 166260.82
## 2015  21536.29  22465.14  25257.55  39269.01  72240.30  82782.06 167201.80
## 2016  22477.27  23406.12  26198.53                                        
##            Aug       Sep       Oct       Nov       Dec
## 2014 157586.43 121528.52  95603.71  42696.74  23377.49
## 2015 158527.41 122469.50  96544.68  43637.72  24318.47
## 2016

次に、放映後ケースについて観光客数を同じ期間について求める。

houeigo.nikka <- c(window(nikka,start=c(2014,4),end=c(2015,3)),fr1.nikka$mean)　# 元のデータから2014年度を切りとり、2015年度予測とつなげる
houeigo.nikka <- ts(houeigo.nikka,freq=12,start=c(2014,4)) # 上でつなげた系列に日付を改めてつける
houeigo.nikka
##           Jan      Feb      Mar      Apr      May      Jun      Jul
## 2014                            18006.00 36305.00 43783.00 58329.00
## 2015 24053.00 41541.00 57506.00 31714.62 41376.21 43783.00 58329.00
## 2016 24053.00 41541.00 57506.00                                    
##           Aug      Sep      Oct      Nov      Dec
## 2014 60355.00 62552.00 74306.00 55082.00 23069.00
## 2015 60355.00 62552.00 74306.00 55082.00 23069.00
## 2016
houeigo.others <- c(window(others,start=c(2014,4),end=c(2015,3)),fr1.others$mean)
houeigo.others<-ts(houeigo.others,freq=12,start=c(2014,4))
houeigo.others
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2014                                25820.00  48997.00  56048.00 125609.00
## 2015  22639.00  25376.00  31192.00  32367.22  53214.50  57818.96 129390.97
## 2016  22579.85  23524.38  27748.73                                        
##            Aug       Sep       Oct       Nov       Dec
## 2014 113332.00 104209.00  72862.00  45846.00  22060.00
## 2015 118695.45 101274.34  72674.15  41555.85  22333.51
## 2016
houeigo.total <- houeigo.nikka+houeigo.others
houeigo.total
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2014                                43826.00  85302.00  99831.00 183938.00
## 2015  46692.00  66917.00  88698.00  64081.84  94590.71 101601.96 187719.97
## 2016  46632.85  65065.38  85254.73                                        
##            Aug       Sep       Oct       Nov       Dec
## 2014 173687.00 166761.00 147168.00 100928.00  45129.00
## 2015 179050.45 163826.34 146980.15  96637.85  45402.51
## 2016

朝ドラ放映の効果は放映後ケースから標準ケースを差し引くことで求められる。

kouka.nikka <- houeigo.nikka - hyoujun.nikka
kouka.others <- houeigo.others - hyoujun.others
kouka <- houeigo.total-hyoujun.total
kouka.nikka
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2014                                2400.608 10926.395 10586.028 17011.633
## 2015 14250.324 28460.156 42931.596 15168.249 15056.625  9645.049 16070.654
## 2016 13309.345 27519.178 41990.618                                        
##            Aug       Sep       Oct       Nov       Dec
## 2014 20328.760 26066.313 40256.584 35303.838 11786.767
## 2015 19387.782 25125.334 39315.605 34362.860 10845.788
## 2016
kouka.others
##             Jan        Feb        Mar        Apr        May        Jun
## 2014                                   3097.4317  3076.3205  7403.9083
## 2015 10905.3846 15991.6996 20508.8533  9644.5837  7293.7870  9174.8479
## 2016 10846.2295 14140.0748 17065.5834                                 
##             Jul        Aug        Sep        Oct        Nov        Dec
## 2014   665.5517 -4228.1883 19166.1674 11307.7114 22927.4181  9964.7395
## 2015  4447.5089  1135.2574 16231.5021 11119.8601 18637.2641 10238.2521
## 2016
kouka
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2014                                5498.039 14002.715 17989.936 17677.184
## 2015 25155.708 44451.856 63440.450 24812.833 22350.412 18819.897 20518.163
## 2016 24155.575 41659.252 59056.202                                        
##            Aug       Sep       Oct       Nov       Dec
## 2014 16100.572 45232.480 51564.295 58231.257 21751.506
## 2015 20523.039 41356.836 50435.465 53000.124 21084.040
## 2016

放映効果はこのようにまず月次ベースで求められた。

年度ベースの放映効果

ニッカについては以下のようになる。

houeigo.nikka.14 <- sum(window(houeigo.nikka,start=c(2014,4),end=c(2015,3)))
houeigo.nikka.15 <- sum(window(houeigo.nikka,start=c(2015,4),end=c(2016,3)))
hyoujun.nikka.14 <- sum(window(hyoujun.nikka,start=c(2014,4),end=c(2015,3)))
hyoujun.nikka.15 <- sum(window(hyoujun.nikka,start=c(2015,4),end=c(2016,3)))
kouka.nikka.14<-sum(window(kouka.nikka,start=c(2014,4),end=c(2015,3)))
kouka.nikka.15<-sum(window(kouka.nikka,start=c(2015,4),end=c(2016,3)))
houeigo.nikka.14
## [1] 554887
houeigo.nikka.15
## [1] 573666.8
hyoujun.nikka.14
## [1] 294578
hyoujun.nikka.15
## [1] 305869.7
kouka.nikka.14
## [1] 260309
kouka.nikka.15
## [1] 267797.1

ニッカ以外の施設については以下のように計算する。

houeigo.others.14 <- sum(window(houeigo.others,start=c(2014,4),end=c(2015,3)))
houeigo.others.15 <- sum(window(houeigo.others,start=c(2015,4),end=c(2016,3)))
hyoujun.others.14 <- sum(window(hyoujun.others,start=c(2014,4),end=c(2015,3)))
hyoujun.others.15 <- sum(window(hyoujun.others,start=c(2015,4),end=c(2016,3)))
kouka.others.14<-sum(window(kouka.others,start=c(2014,4),end=c(2015,3)))
kouka.others.15<-sum(window(kouka.others,start=c(2015,4),end=c(2016,3)))
houeigo.others.14
## [1] 693990
houeigo.others.15
## [1] 703177.9
hyoujun.others.14
## [1] 573203
hyoujun.others.15
## [1] 573203.1
kouka.others.14
## [1] 120787
kouka.others.15
## [1] 129974.8

最後に観光客数合計について計算する。

houeigo.total.14 <- sum(window(houeigo.total,start=c(2014,4),end=c(2015,3)))
houeigo.total.15 <- sum(window(houeigo.total,start=c(2015,4),end=c(2016,3)))
hyoujun.total.14 <- sum(window(hyoujun.total,start=c(2014,4),end=c(2015,3)))
hyoujun.total.15 <- sum(window(hyoujun.total,start=c(2015,4),end=c(2016,3)))
kouka.14<-sum(window(kouka,start=c(2014,4),end=c(2015,3)))
kouka.15<-sum(window(kouka,start=c(2015,4),end=c(2016,3)))
houeigo.total.14
## [1] 1248877
houeigo.total.15
## [1] 1276845
hyoujun.total.14
## [1] 867781
hyoujun.total.15
## [1] 879072.9
kouka.14
## [1] 381096
kouka.15
## [1] 397771.8

以上のように、余市町観光客入込数ベースで2014年度には約38.1万人の上振れがあったと推計され、これは朝ドラ「マッサン」放映の効果と考えられる。この効果は2015年度にも39.8万人程度の当初比上振れとなって現れるものと予想できる。

　【参考】

以下は人数ベースの予測を金額ベースに換算するステップである。詳細は省くので要点だけとする。旅行期間中の平均支出額は北海道観光統計、ニッカ蒸溜所内で実施したアンケート調査に基づく。

金額評価をすると以下のようになる。 まず北海道全体について

round(kouka.14*30040,1)
## [1] 11448123792
round(kouka.15*30040,1)
## [1] 11949066070

次に余市町ベースで金額換算をする。

round(kouka.14*7875,1)
## [1] 3001130988
round(kouka.15*7875,1)
## [1] 3132453239