BLOG

 

235日のセブIT留学  成長日記「時系列分析_④」 ( 18/235 days )

26 7月 2018, Posted by keisuke in IT留学, プログラミング, 統計

こんばんわ!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): #標本自己共分散ーーーパラメーターφ1

   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に実装

・次数の打ち込み

Post a comment