こんばんわ！Keisukeです！！
今日は, AR(p=1)モデルを用いてアイスクリームの消費量を計算したいと思います！！
Excelでも計算することはできますが、途方もなく時間がかかってしまいます！
こういう時は, Pythonですね！
コード
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
import csv
from time import time
import numpy as np
import matplotlib.pyplot as plt
def scatter_plot(x,y,x1,y1): #グラフ表示
plt.plot(x,y,marker=’o’)
plt.plot(x1,y1,marker=’x’)
plt.xlabel(‘day’)
plt.ylabel(‘temp’)
plt.show()
def read_csv(filename): #CSV読み込み
X = []
Y = []
with open(filename) as f:
reader = csv.reader(f)
next(reader)
for row in reader:
X.append(float(row[0]))
Y.append(float(row[1]))
return X, Y
def cal_ave(data): #平均値取得
s = sum(data)
N = len(data)
ave = s/N
return ave
def cal_dif(data): #階差
dif = []
for i in range(len(data)):
dif.append(data[i]-data[(i-1)])
dif[0] = 0
return dif
def cal_add(data): #積算
add = []
add.append(0)
for i in range(len(data)):
add.append(data[i]+add[(i)])
add.pop(0)
return add
def cal_add_ave(data):#積算平均
addave = []
add = cal_add(data)
for i in range(len(data)):
addave.append(add[i]/(i+1))
return addave
def cal_dev0(data): #偏差(y-ave)^2
dev0 = []
addave = cal_add_ave(data)
for i in range(len(data)):
for n in range(i+1):
dev0.append((data[n]-addave[i])**2)
return dev0
def cal_dev_sum0(data): #合計偏差(y-ave)^2
dev0 = cal_dev0(data)
dev0_sum = []
#print(‘Start’)
for i in range(len(data)):
dev0_sum_step = 0
#print(dev)
for n in range(i+1):
dev0_sum_step = dev0_sum_step + dev0[n]
for i in range(i+1):
dev0.pop(0)
#print(len(dev))
dev0_sum.append(dev0_sum_step/(i+1))
#print(‘dev_sum_step’)
#print(dev_sum)
return dev0_sum
def cal_dev1(data): #偏差(y – ave)(y-1 – ave)
dev1 = []
addave = cal_add_ave(data)
data0 = []
data0.append(0)
for l in data:
data0.append(l)
#print(data)
#print(data0)
#print(addave)
for i in range(len(data)):
#print(‘start’)
for n in range(i+1):
if i>0 and n>0:
ai = (data[n]-addave[i])*(data0[n]-addave[i])
dev1.append(ai)
else:
dev1.append(0.0)
#print(data[n]-addave[i])
#print(data0[n]-addave[i])
#print(ai)
#print(dev1)
return dev1
def cal_dev_sum1(data): #合計偏差(y – ave)(y-1 – ave)
dev1 = cal_dev1(data)
dev1_sum = []
#print(‘Start’)
for i in range(len(data)):
dev1_sum_step = 0
#print(dev1)
for n in range(i+1):
dev1_sum_step = dev1_sum_step + dev1[n]
for i in range(i+1):
dev1.pop(0)
#print(len(dev1))
dev1_sum.append(dev1_sum_step/(i+1))
#print(‘dev1_sum_step’)
#print(dev1_sum)
return dev1_sum
def cal_cor(data): #標本自己共分散ーーーパラメーターφ１
cor = []
dev_sum0 = cal_dev_sum0(data)
dev_sum1 = cal_dev_sum1(data)
for i in range(len(data)):
if dev_sum0[i]>0:
cor.append(dev_sum1[i]/dev_sum0[i])
else:
cor.append(0)
return cor
def cal_pred(data): #予測値
pred = []
#pred.append(0)
addave = cal_add_ave(data)
cor = cal_cor(data)
#print(‘addave’)
#print(addave)
#print(‘cor’)
#print(cor)
data0 = []
data0.append(0)
for l in data:
data0.append(l)
for i in range(len(data)):
y = cor[i]*(data0[i]-addave[i])+addave[i]
pred.append(y)
print(cor[10]*(data0[8]-addave[9])+addave[9])
return pred
if __name__ == ‘__main__’:
tm = time() # Timer Start
data0, data1 = read_csv(‘ice.csv’)
N = len(data0)
sum0 = sum(data0)
sum1 = sum(data1)
ave0 = cal_ave(data0)
ave1 = cal_ave(data1)
data2 = cal_dif(data1)
data3 = cal_add(data1)
data4 = cal_add_ave(data1)
data5 = cal_dev0(data1)
data6 = cal_dev_sum0(data1)
data7 = cal_dev1(data1)
data8 = cal_dev_sum1(data1)
data9 = cal_cor(data1)
data10 = cal_pred(data1)
scatter_plot(data0, data1, data0, data10)
print(“Runtime : %.3f [sec]”%(time()-tm)) # Timer Stop & Disp
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
実行結果
うんん～、一月のデータだけうまく行っていない気がする…そのうち修正します（笑）
時系列分析の各モデルの計算はデータ数が増えれば増えるほど多量な計算コストがかかってしまいます！Yule-Waker方程式を見てもらうと一目瞭然だと思います！
こういう計算は, プログラミングでやっちゃうのが効率いいですね！
随時, 理解が追い付いたらMAモデル, ARMAモデル, ARIMAモデル, SARIMAモデルを記述してきます！
ARモデルも次数打ち込んで対応できるようにしたいな….
[今日の達成]
・AR(p=1)モデルをPythonに実装
[今日の未消化]
・バグ？がありそう
・他のモデルをPythonに実装
・次数の打ち込み