【Python】統計学的モデル　時系列分析の実装

#必要なライブラリをインポートする　#対数変換を行う
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import numpy as np

#データを読み込む
# データの読み込み(StatsModelsの検証データを用います)
macrodata = sm.datasets.macrodata.load_pandas().data
macrodata.index = pd.Index(sm.tsa.datetools.dates_from_range("1959q1","2009q3"))

# 読み込んだデータからアメリカの実質GDPの対数変換前の値を表示します。
print(macrodata.realgdp.head())

#macrodataを対数変換する
macrodata_realgdp_log = np.log(macrodata.realgdp)

print(macrodata_realgdp_log.head())

#必要なライブラリをインポートし、階差系列をとってみる
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
import numpy as np

#statsmodelの検証データを読み込む
co2_tsdata = sm.datasets.co2.load_pandas().data

#欠損値の処理
co2_tsdata2 = co2_tsdata.fillna(method="ffill")

#データの階差をとる
co2_tsdata2_diff = co2_tsdata.diff()

#階差前と階差後のグラフを可視化、比較する
plt.subplot(2,1,1)
plt.title("Mauna Loa Weekly Atmospheric CO2 Data")
plt.xlabel("date")
plt.ylabel("CO2 Concentration ppmv")
plt.plot(co2_tsdata2)

plt.subplot(2,1,2)
plt.title("Mauna Loa Weekly Atmospheric CO2 Data DIFF")
plt.xlabel("date")
plt.ylabel("CO2 Concentration ppmv DIFF")
plt.plot(co2_tsdata2_diff)

plt.subplots_adjust(wspace=0, hspace=1.0)

plt.show()

#季節調整済系列、seasonal_decomposeでトレンド・季節変動・不規則変動に分けることができる
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from numpy import nan as na
import numpy as np

# データの読み込み(StatsModelsの検証データを用います)
co2_tsdata = sm.datasets.co2.load_pandas().data
co2_tsdata2 = co2_tsdata.fillna(method="ffill")
# 季節調整を行い原系列をトレンド、季節変動、残差に分けて出力します
fig = sm.tsa.seasonal_decompose(co2_tsdata2, period=52).plot()
plt.show()

#移動平均での変換
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from datetime import datetime
import numpy as np

co2_tsdata = sm.datasets.co2.load_pandas().data

# 欠損値の処理
co2_tsdata2 = co2_tsdata.fillna(method="ffill")

# 原系列のグラフ
plt.subplot(6, 1, 1)
plt.xlabel("date")
plt.ylabel("co2")
plt.plot(co2_tsdata2)

#移動平均を求める
co2_moving_avg = co2_tsdata2.rolling(window=51).mean()
#移動平均のグラフを描画する
plt.subplot(6,1,3)
plt.xlabel("date")
plt.ylabel("co2")
plt.plot(co2_moving_avg)

#原系列-移動平均のグラフを描画する
plt.subplot(6, 1, 5)
plt.xlabel("date")
plt.ylabel("co2")
mov_diff_co2_tsdata = co2_tsdata2-co2_moving_avg 
plt.plot(mov_diff_co2_tsdata)
plt.show()

#必要なライブラリをインポート、データを読み込む
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm

url = 'https://www.salesanalytics.co.jp/l6p7'
df_wine = pd.read_csv(url)

#1980年1月から1994年9月までの月次のインデックスを作成し、月毎の頻度を指定。
#生成した日付インデックスを、dw_wineのインデックスにせていし、元のMonthlyの列を削除する
#これらの処理で、月次売上データを扱うデータフレームを準備している
index = pd.date_range("Jan-80","Sep-94",freq="M")
df_wine.index=index
del df_wine["Monthly"]

# 最適なパラメータを算出するため関数の定義
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)]

# 最適なパラメーターを算出
selectparameter(df_wine,12)

#自己相関、偏自己相関を可視化する
fig=plt.figure(figsize=(12, 8))

# 自己相関係数のグラフを出力します
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df_wine, lags=30, ax=ax1)

# 偏自己相関係数のグラフを出力します
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df_wine, lags=30, ax=ax2)
plt.show()

#SARIMAモデルを構築する
SARIMA_df_wine = sm.tsa.statespace.SARIMAX(df_wine, order=(0,0,0), seasonal_order=(0,1,0,12)).fit()

#構築したSARIMAモデルのBICを出力します
print(SARIMA_df_wine.bic)

#predに予測データを代入する
pred = SARIMA_df_wine.predict("1994-7-31", "1997-12-31")

#SARIMAで予測したモデルと元のデータの可視化
plt.plot(pred, c="r")
plt.plot(df_wine)
plt.show()

#必要なライブラリをインポート、データを読み込む
import warnings
import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm

url = 'https://www.salesanalytics.co.jp/l6p7'
df_wine = pd.read_csv(url)

#月次売上データを扱うデータフレームを準備
index = pd.date_range("Jan-80","Sep-94",freq="M")
df_wine.index=index
del df_wine["Monthly"]

# 最適なパラメータを算出するため関数の定義
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)]

# 最適なパラメーターを算出
selectparameter(df_wine,12)

#自己相関、偏自己相関を可視化する
fig=plt.figure(figsize=(12, 8))

# 自己相関係数のグラフを出力します
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(df_wine, lags=30, ax=ax1)

# 偏自己相関係数のグラフを出力します
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(df_wine, lags=30, ax=ax2)
plt.show()

#SARIMAモデルを構築する
SARIMA_df_wine = sm.tsa.statespace.SARIMAX(df_wine, order=(0,0,0), seasonal_order=(0,1,0,12)).fit()

#構築したSARIMAモデルのBICを出力します
print(SARIMA_df_wine.bic)

#predに予測データを代入する
pred = SARIMA_df_wine.predict("1994-7-31", "1997-12-31")

#SARIMAで予測したモデルと元のデータの可視化
plt.plot(pred, c="r")
plt.plot(df_wine)
plt.show()