初めまして。Qiita初投稿・まだまだ勉強中のばりきゃりを目指しつつ適度に働きたい（@moku0517)です。
学習のアウトプット発信の一環として本記事を投稿します。まだまだ勉強中のみですので、暖かくご指摘・コメントいただければ幸いです。
本分析の目的は、旅行会社で働く私が、避暑地ってだんだん避暑地じゃなくなるのか？そうなったら観光客は来てくれるのか？それとも新たな魅力を発信しないとだめなのか？という疑問を解消するためにおこないます。
旅行会社ですが地域脱炭素普及のための事業をしていることもあり、このような分析をしてみたいと思いました。
- 過去18年の軽井沢の観光客数の変遷と、温室効果ガス排出量の合計をそれぞれデータにする
- 相関係数を確認
- 今後の観光客数を予測
①温室効果ガスの排出量と軽井沢の観光客数との関連を調べる
日本の温室効果ガスのデータを読み込みます。
from google.colab import files
uploaded = files.upload()
データが読み込めたら、そのExcelの中の、Summaryというシートを取り出します。
import pandas as pd
xlsx = pd.ExcelFile("./L5-7gas_2023_gioweb_ver1.0.xlsx")
Summary = pd.read_excel(xlsx,"1.Summary")
これを一回出してみると。。
列も行も合っていないため、データがうまく取り出せませんでした。
そこで、行を揃え、列は列番号を指定して取り出したいデータを取り出せるようにしました。
import pandas as pd
xlsx = pd.ExcelFile("./L5-7gas_2023_gioweb_ver1.0.xlsx")
Summary = pd.read_excel(xlsx,"1.Summary",skiprows=3,usecols=list(range(26,58)))
こうすると、下記のように正しいデータが取り出せています。
今回は、10行目までの数値を取り出したかったので、10行目までを指定すると下記のようになったので、下記のように変数を指定します。
GHG_Summary=Summary.loc[10]
今回は、軽井沢の観光客数との相関を調べるのですが、該当するデータが2002年からなので、データの個数を合わせるため、下記の画像の通り出力するための計算式を、指定します。
GHG_Summary_2002=GHG_Summary.loc[2002:]
これで、日本の温室効果ガスの排出量合計の2002年から2021年までの変遷が分かるようになりました。
- ２．軽井沢の観光客数のデータを取り出す
今回参考にする、軽井沢の観光客数の元データは、長野県が出している各地域の観光客数の推移のデータでした。
画像の通り、データのしにくいので、手作業で該当するデータ（2022年から、2022年までのデータ）を取り出し、下記の通りデータフレームを作成します。
num = [8051,7635,7772,7737,7820,7915,7691,7630,7759,7701,7796,7946,8277,8403,8458,8530,8707,8423,5144,5484,7129]
travel = pd.DataFrame(num, index = range(2002,2023), columns=["観光客数"])
最初に作成した日本の温室効果ガスの排出量のデータ量と一致させるため、2021年までのデータを使います。
-1にすると、もともとあったデータのうち、最後の行は無くなります。
travel[:-1]
これを改めてtravelとしておきます。
travel=travel[:-1]
ここまでで、軽井沢の観光客数のデータを取り出せました、
次に、2つのデータを組み合わせます。
- ３．日本の温室効果ガスの排出量と軽井沢の観光客数のデータを組み合わせる
組み合わせる時は下記のコードとします。
travel["GHG_Summary"]=GHG_Summary_2002
これを実行すると。。
ついに一つのデータとすることができました！
上記まで来たらあとは相関係数を求めるだけです！
相関係数は下記の通り求めます。
travel.corr()
実行すると。。。ついに出すことができました！
一見、正の相関があるようですが、どうでしょうか。。？
一見、温室効果ガスの排出量が増えると観光客数も増えているように見えます。
この時点で私が求めたかった、温室効果ガスの排出量が増えると気温も上がるから、避暑地の観光客数は少しずつ減るのではないか？という推測は現時点では当てはまらないことがわかりました。
ただ、一方で、これは、疑似相関ではないかとも考えられます。
観光客数が増えることと、温室効果ガスの排出量が増えることは、確かに、人の移動が増えることで車や店舗等からの排出量が増えるのでつながりがありそうですが、例えば、観光客のボーナスが年々上がっているであったり、軽井沢が新しいホテルを建設したから増えているなどの、他の理由も考えられそうです。
データの使い方は難しいのですね。。
もしかしたら、私が調べたかったことは、扱うデータが違えばまた違う答えが出るのかもしれません。
②観光客数を予測する（モデル構築（SARIMAモデル））
ここからは、少々視点を変えて、今後の軽井沢の観光客数を予測するモデルを構築したいと思います。
今回は、SARIMAモデルで構築します。
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
%matplotlib inline
今回、データは、長野県の統計資料を使用しました。
長野県データ
予測には、上述の通り、SARIMAモデルを使いたかったので、各年の1月から12月までのデータを取り出す必要があり、この資料からだとかなりデータがとりづらかったのですが、PDFから必要なデータを取り出し、CSVにしました。（平成15年から令和4年までの観光客数）
CSVを読み込むと、下記の表のようになります。
df = pd.read_csv("./軽井沢観光者数月別まとめ.csv")
df.head()
- ２．データの形をモデル構築用に変形する
上記ではまだSARIMAモデルには使えないので、データの形を変えたいと思います。
まずは、日付のカラムを取得します。
今回はデータをまずリスト型で抽出したいので、格納する場所を用意しておきます。
y_m_list = list()
data_list = list()
for row in df.iterrows():
_row = row[1].values
_cols = row[1].index.values
続いて、今回元データでは、年月（日）が和暦になっていたので、それを西暦に直しておくことが後々大事になってきます。
import unicodedata
import re
import datetime
# 各年号の元年を定義
eraDict = {
"明治": 1868,
"大正": 1912,
"昭和": 1926,
"平成": 1989,
"令和": 2019,
}
def japanese_calendar_converter(text):
# 正規化
normalized_text = unicodedata.normalize("NFKC", text)
# 年月日を抽出
pattern = r"(?P<era>{eraList})(?P<year>[0-9]{{1,2}}|元)年(?P<month>[0-9]{{1,2}})月(?P<day>[0-9]{{1,2}})日".format(eraList="|".join(eraDict.keys()))
date = re.search(pattern, normalized_text)
# 抽出できなかったら終わり
if date is None:
print("Cannot convert to western year")
# 年を変換
for era, startYear in eraDict.items():
if date.group("era") == era:
if date.group("year") == "元":
year = eraDict[era]
else:
year = eraDict[era] + int(date.group("year")) - 1
# date型に変換して返す
return datetime.date(year, int(date.group("month")), int(date.group("day")))
ここまで準備できたら、ついに日付のカラムを取得します。
# 日付のカラムを取得
year = _row[0]
for i, month in enumerate(_cols):
if i == 0: continue
_y_m = japanese_calendar_converter(f"平成{year}年{month}1日")
y_m_list.append(_y_m)
出力するとdatetime.date(2003, 1, 1),・・・にように2022,12,1まで出てくると思います。
これでいったん成功です。
次に、観光客数のデータを取得します。
基本的に上記と同様のやり方です。
後段は、数字に入っている、カンマを取り除いたり、文字列で入っているデータを氷解し元のデータに戻す作業をします。
for i, data in enumerate(_row):
if i == 0: continue
_data = data.replace(",", "")
_data = eval(_data)
data_list.append(_data)
出力すると2498から4105までのリストが出てきます。
ここまで出来たら後はデータとして成型するだけ。
traveler = pd.DataFrame({"Month": y_m_list, "sales": data_list})
データの初めと終わりを出力してみます。
欲しかったデータを取得できました。
データの前処理ができたら、モデルを構築します。
私がデータ分析を学んだAidemyさんで学習した時のプログラムを書いていきます。
# travelerのインデックスにindexを代入
traveler.index = traveler["Month"]
# travelerの"Month"カラムを削除
del traveler["Month"]
# orderの最適化関数
def selectparameter(DATA, s):
p = d = q = range(0, 2)
pdq = list(itertools.product(p, d, q))
seasonal_pdq = [(x[0], x[1], x[2], s) for x in list(itertools.product(p, d, q))]
parameters = []
BICs = np.array([])
for param in pdq:
for param_seasonal in seasonal_pdq:
try:
mod = sm.tsa.statespace.SARIMAX(DATA,
order=param,
seasonal_order=param_seasonal)
results = mod.fit()
parameters.append([param, param_seasonal, results.bic])
BICs = np.append(BICs, results.bic)
except:
continue
return parameters[np.argmin(BICs)]
理解しきれていない部分もあり恐縮ですが。。。
上記は情報量規準 (今回の場合は BIC(ベイズ情報量基準) )によってどの値が最も適切なのか調べるプログラムです。
各パラメーターがそれぞれ、0,1の場合についてのSARIMAモデルのBICを計算し、最もBICが小さくなった場合を表示しています。
最適化ができたら、ついにモデルの構築です。周期は月ごとのデータであることも考慮してs=12とします。
best_params = selectparameter(traveler, 12)
SARIMA_traveler = sm.tsa.statespace.SARIMAX(traveler,order=(0,0,0),seasonal_order=(0,1,1,12)).fit()
いよいよ最後のプログラムです。
今までのモデルから、今回は、2022年12月1日～2024年12月末までの観光客数を予測します。
pred = SARIMA_traveler.predict("2022-12-01", "2024-12-29")
グラフ化します。予測部分は赤で記載します。
plt.plot(traveler)
plt.plot(pred,color="r")
すると、下記のような画像が出力されました。
上記のようになりました。
実際は、コロナ禍のデータも入ってしまっているので、今後はもう少し増えるのかなとも思いつつ、一定程度、1年に1回ピークが来ていることがわかるデータとなりました。
今後は、気温との変化も注視しつつ、観光客の推移をより分析できるようにしていきたいと思います。
拙い文章、分析にお付き合いいただきまして、ありがとうございました。
超初心者としては、頑張ったほうかと思っていますが、使用する元データが適切でなかったり、成果物までの時間がなく簡単にしてしまったりと反省すべき点は多くあります。
今回はここまでで終了しますが、上記に記載している通り、もう少し自分の知りたいことに近づけられるような分析ができるといいなと思っておりますので、引き続き勉強していきたいです。
長々とお付き合いいただきありがとうございました！