はじめに

Rのforecastパッケージのetsを使用して、AirPassengersの直近3年分から指数平滑法モデルを作成。
平滑係数を変化させてAICのグラフを描いてみると、微分係数がゼロでないものがあり、まだ改善の余地がある。
etsのETS(A)モデルをtensorflowで実行してみた。
tensorflowで最適化したら、etsよりも二乗誤差が小さい解が得られた。但し、平滑係数はすべてゼロ!?
制約条件の無いパラメータの微分係数はほぼゼロとなっているため、少なくとも局所解となっているはず。
forecastパッケージのソースを見ても、ETS(A)の場合、二乗誤差の対数を最小化しているように見える。

Rスクリプト

library(forecast)

n=36
history=AirPassengers[(144-n+1):144]

m1=mean(history)
m2=sd(history)
y1=(history-m1)/m2
ds1 = ts(y1,start=c(1958,1),fr=12)
fmodel1 = ets(ds1)

etsによるモデル

ETS(A,N,A) 

Call:
 ets(y = ds1) 

  Smoothing parameters:
    alpha = 0.649 
    gamma = 0.0216 

  Initial states:
    l = -0.4015 
    s=-0.7145 -1.1136 -0.3391 0.3096 1.6001 1.6217
           0.7145 -0.0071 -0.2135 -0.3195 -0.8903 -0.6485

  sigma:  0.1515

     AIC     AICc      BIC 
23.14316 47.14316 46.89595

指数平滑法ETS(A,A,A)の更新式

\begin{align}
l_t&=\alpha(y_t-s_{t-m})+(1-\alpha)(l_{t-1}+b_{t-1})\\
b_t&=\beta^*(l_t-l_{t-1})+(1-\beta^*)b_{t-1}\\
s_t&=\gamma(y_t-l_{t-1}-b_{t-1})+(1-\gamma)s_{t-m}\\
p_t&=l_{t-1}+b_{t-1}+s_{t-m},\quad e_t=y_t-p_t\\
l_t&=l_{t-1}+b_{t-1}+\alpha e_t\\
b_t&=b_{t-1}+\alpha\beta^*e_t\\
s_t&=s_{t-m}+\gamma e_t\\
AIC&=L^*+2q,\quad L^*=n\log\left(\sum e_t^2\right)
\end{align}

Pythonスクリプト

import numpy as np
import pandas as pd
import tensorflow as tf

# Exponential smoothing ETS(A,A,A)
def ets(history, params):
  alpha = params[0]
  beta  = params[1]
  gamma = params[2]
  level = params[3]
  trend = params[4]
  seasonal = params[5:16]
  seasonal.append(-sum(seasonal))
  predict = []
  error = []
  n = len(history)
  for i in range(n):
    s = i % 12
    pbase = level + trend
    predict.append(pbase + seasonal[s])
    error.append(history[i] - predict[i])
    level = pbase + alpha * error[i]
    trend = trend + beta * error[i]
    seasonal[s] = seasonal[s] + gamma * error[i]
  return error

# Exponential smoothing ETS(A,A,A)
def tfEts(history, n, params):
  alpha = params[0]
  beta  = params[1]
  gamma = params[2]
  level = params[3]
  trend = params[4]
  seasonal = tf.split(params[5:16], 11)
  s11 = params[5:16]
  s12 = tf.negative(tf.reduce_sum(s11))
  seasonal.extend([s12])
  predict = []
  error = []
  for i in range(n):
    s = i % 12
    pbase = level + trend
    predict.append(pbase + seasonal[s])
    error.append(history[i] - predict[i])
    level = pbase + alpha * error[i]
    trend = trend + beta * error[i]
    seasonal[s] = seasonal[s] + gamma * error[i]
  err = tf.stack(error)
  return err

if __name__ == '__main__':
  AirPassengers = [
    340,318,362,348,363,435,491,505,404,359,310,337,
    360,342,406,396,420,472,548,559,463,407,362,405,
    417,391,419,461,472,535,622,606,508,461,390,432]
  history = np.array(AirPassengers)
  m1 = history.mean()
  m2 = history.std(ddof=1)
  y1 = (history - m1) / m2
  N  = y1.size

  # etsにより得られたパラメータ
  params = [0.6489863, 0, 0.0216017]
  params.extend([-0.4014818, 0])
  params.extend([
    -0.64847373, -0.89031436, -0.31948061, -0.21345777,
    -0.00713472,  0.71451547,  1.62174834,  1.60011895,
     0.30963273, -0.33910409, -1.11356078])

  error = ets(y1, params)
  sse = sum([e * e for e in error])
  print("sse=%f" % (sse / N))
  print("AIC=%f" % (N * np.log(sse) + 2 * 15))

  # 回帰直線を求めてパラメータを初期化
  x = np.arange(y1.size)
  x = x - x.mean()
  b = sum(x * y1) / sum(x * x)
  l = y1 - b * x

  params = [0.0, 0.0, 0.0]
  params.extend([(x[0] - 1) * b, b])
  params.extend([l[i] for i in range(11)])

  yt = tf.placeholder(tf.float32, shape=[N])
  pt = tf.Variable(params[0:16], dtype=tf.float32)
  err = tfEts(yt, N, pt)
  loss = tf.reduce_mean(tf.square(err))
  grad = tf.gradients(loss, pt)
  train_step = tf.train.AdamOptimizer(1e-3).minimize(loss)
  pt_ = tf.clip_by_value(pt, [0]*3+[-100]*13, [1]*3+[100]*13)
  clip_step = tf.assign(pt, pt_)
  result = []
  CYCLE = 1000
  MD = 100
  with tf.Session() as sess:
    sess.run(tf.global_variables_initializer())
    for i in range(CYCLE):
      train_step.run(feed_dict={yt: y1})
      sess.run(clip_step)
      if (i+1) % MD == 0:
        a,b,c = sess.run([loss,pt,grad], feed_dict={yt: y1})
        result.append(np.concatenate([[a],b,c[0]]))
        print("i=%d,loss=%f" % (i, a))

  names = ['loss']
  names.extend(["alpha","beta","gamma","level","trend"])
  names.extend(['s'+str(i) for i in range(11)])
  names.extend(["alpha","beta","gamma","level","trend"])
  names.extend(['d'+str(i) for i in range(11)])
  df = pd.DataFrame(result, columns=names)
  df.to_csv("result.txt", sep="\t", index=False, encoding="shift_jis")

  print("%s, %f" % (names[0], a))
  for i in range(16):
    print("%s, %f, %f" % (names[i+1], b[i], c[0][i]))

  params = list(b)
  error = ets(y1, params)
  sse = sum([e * e for e in error])
  print("sse=%f" % (sse / N))

実行結果

sse=0.022960
AIC=23.143163
i=99,loss=0.018410
i=199,loss=0.015501
i=299,loss=0.014946
i=399,loss=0.014882
i=499,loss=0.014878
i=599,loss=0.014877
i=699,loss=0.014877
i=799,loss=0.014877
i=899,loss=0.014877
i=999,loss=0.014877
loss, 0.014877
alpha, 0.000000, 0.014877
beta, 0.000000, 0.127179
gamma, 0.000000, 0.014877
level, -0.924725, 0.000001
trend, 0.049985, 0.000021
s0, -0.433102, 0.000000
s1, -0.760412, 0.000000
s2, -0.238940, -0.000000
s3, -0.213290, -0.000000
s4, -0.053180, -0.000000
s5, 0.682590, -0.000000
s6, 1.552820, -0.000000
s7, 1.540654, -0.000000
s8, 0.251107, -0.000000
s9, -0.420759, -0.000000
s10, -1.164059, -0.000000
sse=0.014877