\( \DeclareMathOperator{\abs}{abs} \newcommand{\ensuremath}[1]{\mbox{$#1$}} \)
Enders 1995 Applied Econometric Time Series の第1章では
2次差分方程式（定係数、同次）で特性方程式が複素根を持つ場合,
一般解は、最終的には三角関数で表現できることが示されている。
これは三角関数の加法定理や倍角公式、ド・モアブル定理などを使って式を
変形すると共役複素数の解を導き、その和を求めると虚数部分が消去されて
cos()コサインの実関数が導出される。エンダースのテキストではこの導出過程が示されているが、
数学ソフトでは下記に示すように複素数で表示されておりcos()の解を表示してくれない。
AIに一発でcos()の型の解を導く方法はないか聞いたが、今のところは、
人間の手作業で式を整理してコサインの式を導かざるを
得ないようである。
Endersの第1章の補論1の式(A1.11)を第1章のworksheet1.2
のcase1 の差分方程式を例として数学ソフトMAXIMAを使って、コサインの関数形に導く
計算実験してみた。
case 1 では y(t)=1.6y(t-1)-0.9y(t-2) の2次の線形差分方程式を解いている。
テキストでは初期値は与えられていないが、ここではy(0)=1, y(1)=-0.8 と仮定する。
差分方程式については和書では参考文献、小山 1995 数学教室で解説されている。
線形の2次差分方程式（定係数、同次）の一般解は特解（基本解と呼ばれることもある）の
線形結合で与えられる。y(t)+a*y(t-1)+b*y(t-2)=0のように方程式が恒等的にゼロとなる場合は
同次と呼ばれる。y(t)+a*y(t-1)+b*y(t-2)=R(t)のような非同次の場合には同次の場合の一般解に
未定係数法や演算子法で求めた特殊解の線形結合を一般解と呼んでいる。用語が紛らわしいので
同次の場合の一般解を同次解と呼び、それに非同次の場合の特殊解を加えたものが非同次差分方程式
の一般解と呼ばれている。参考文献の数学書では非同次方程式の右辺をゼロとおいて、これを
余方程式(complementary equation)とし、余方程式の一般解を余関数と呼び、非同次線形差分方程式の
一般解を特殊解と余関数の和と説明されている。この余関数は同次解と同じ意味となる。
Endersの第１章での差分方程式の解法は
４段階で説明されている。第１に同時方程式を設定し、同時解を得る。第２に特殊解を１つ見つける。
第３に一般解を特殊解を全ての同時解の線形結合の和とする。第４に一般解に初期条件の制約を
課して任意の係数を特定する。いずれも非同次の差分方程式の解法を語っており
エンダースが示す４段階の作業は y(t)+a*y(t-1)+b*y(t-2)=R(t)のような非同次の差分方程式
の解法を示している。
以下の計算例ではy(t)-1.6y(t-1)+0.9y(t-2)=0 という２次の同次差分方程式を使っている。
| --> | /* 邦訳書 第1章 補論(A1.11 ) 1995原書 APPENDIX 1 (A1.11) */; |
差分方程式を解くためのパッケージをロードしておく。
| (%i1) |
load("solve_rec"); |
\[\]\[\tag{%o1} "C:/maxima-5.49.0/share/maxima/5.49.0/share/solve\_ rec/solve\_ rec.mac"\]
| (%i2) | ratprint: false$ |
Endersのテキストでは初期値は与えられていないがy(0)=1,y(1)=-0.8と仮定して話しを進める。
| (%i3) | sol:solve_rec(y[t]=1.6·y[t−1]−0.9·y[t−2],y[t],y[0]=1,y[1]=−0.8) ; |
\[\]\[\tag{sol} {y_t}\mathop{=}\frac{{{\left( \sqrt{26} \% i\mathop{+}8\right) }^{t}} \left( 8 \sqrt{26} \% i\mathop{+}13\right) }{26 {{10}^{t}}}\mathop{-}\frac{{{\left( 8\mathop{-}\sqrt{26} \% i\right) }^{t}} \left( 8 \sqrt{26} \% i\mathop{-}13\right) }{26 {{10}^{t}}}\]
上記のy(t)式をrectformを使ってxy座標表示に変えてると、教数部分は消えている。これを改めてy(t)の関数で定義する。
| (%i4) | rectform(sol),float; |
\[\]\[\tag{%o4} {y_t}\mathop{=}\mathop{(}0.038461538461538464 \mathop{(}13.0 {{9.486832980505138}^{t}} \cos{\left( 0.5674504752788371 t\right) }\mathop{-}40.792156108742276 {{9.486832980505138}^{t}} \sin{\left( 0.5674504752788371 t\right) }\mathop{)}\mathop{)}/{{10}^{t}}\mathop{-}\]
| (%i5) |
y(t):=(0.038461538461538464·(13.0·9.486832980505138^t·cos(0.5674504752788371·t) −40.792156108742276·9.486832980505138^t·sin(0.5674504752788371·t)))/10^t −(0.038461538461538464·(40.792156108742276·9.486832980505138^t·sin(0.5674504752788371·t) −13.0·9.486832980505138^t·cos(0.5674504752788371·t)))/10^t $ |
y(t)をグラフで表してみると大きな振幅の波動が時間と共に収束している。
| (%i6) |
wxplot2d(y(t),[t,0,50], [style, [lines, 3, 2]]); |
\[\]\[\tag{%t6} \]
\[\]\[\tag{%o6} \]
2次の差分方程式の特性根α1について極形式で表し、絶対値と偏角を求める。
| (%i7) | sol:solve_rec(y[t]=1.6·y[t−1]−0.9·y[t−2],y[t]) ; |
\[\]\[\tag{sol} {y_t}\mathop{=}\frac{{{\left( \sqrt{26} \% i\mathop{+}8\right) }^{t}} {{\ensuremath{\mathrm{\% k}}}_2}}{{{10}^{t}}}\mathop{+}\frac{{{\left( 8\mathop{-}\sqrt{26} \% i\right) }^{t}} {{\ensuremath{\mathrm{\% k}}}_1}}{{{10}^{t}}}\]
| --> | /* α1=a+b i の計算 */; |
| (%i8) | (sqrt(26)·%i+8)/10,float; |
\[\]\[\tag{%o8} 0.1 \left( 5.0990195135927845 \% i\mathop{+}8\right) \]
| (%i9) | z:0.1·5.0990195135927845·%i+0.1·8; |
\[\]\[\tag{z} 0.5099019513592785 \% i\mathop{+}0.8\]
| (%i10) | cabs(z); |
\[\]\[\tag{%o10} 0.9486832980505139\]
α1の絶対値は0.94868と計算できた。
| (%i11) | carg(z); |
\[\]\[\tag{%o11} 0.5674504752788371\]
α1の偏角は0.56745と計算できた。
エンダースのテキストの補論１にあるようにCOSの関数に整理される。
Enders 1995 (1.51),邦訳書 (1.37)の式に絶対値と偏角をあてはめる。
1995原書の式1.51、邦訳書の式1.37 の式が得られる。
| (%i12) | ycos(t):=C1·0.9486832980505139^t·cos(0.567450475278837·t + C2); |
\[\]\[\tag{%o12} \mathop{ycos}(t)\mathop{:=}\ensuremath{\mathrm{C1}} {{0.9486832980505139}^{t}} \cos{\left( 0.567450475278837 t\mathop{+}\ensuremath{\mathrm{C2}}\right) }\]
任意の定数C1,C2を特定さすため初期値y(0)=1,y(1)=-0.8の仮定からC1、C2の 非線型の連立方程式を設定して解いてみる。近似解を得るためにニュートン法で C1、C2を求める。ニュートン法を使うため必要なパッケージをロードしておく。
| (%i15) |
load("mnewton")$ eq1: C1·0.9486832980505139^0·cos(0.5674504752788389·0+C2) = 1; eq2: C1·0.9486832980505139^1·cos(0.5674504752788389·1+C2) = −0.8; |
\[\]\[\tag{eq1} 1.0 \ensuremath{\mathrm{C1}} \cos{\left( \ensuremath{\mathrm{C2}}\right) }\mathop{=}1\]
\[\]\[\tag{eq2} 0.9486832980505139 \ensuremath{\mathrm{C1}} \cos{\left( \ensuremath{\mathrm{C2}}\mathop{+}0.5674504752788389\right) }\mathop{=}\mathop{-}0.8\]
| (%i16) |
/* 数値解の計算 (C1, C2) */ mnewton([eq1, eq2], [C1,C2], [1, 1]); |
\[\]\[\tag{%o16} \left[ \left[ \ensuremath{\mathrm{C1}}\mathop{=}3.293349942862704\mathop{,}\ensuremath{\mathrm{C2}}\mathop{=}1.2622833124845299\right] \right] \]
ycos(t)の式にC1,C2の近似解を代入して、ycos(t)を関数yy(t)定義し直して グラフを描いてみる。
| (%i17) | yy(t):=subst([C1=3.293349942862704,C2=1.2622833124845299],ycos(t)); |
\[\]\[\tag{%o17} \mathop{yy}(t)\mathop{:=}\mathop{subst}\left( \left[ \ensuremath{\mathrm{C1}}\mathop{=}3.293349942862704\mathop{,}\ensuremath{\mathrm{C2}}\mathop{=}1.2622833124845299\right] \mathop{,}\mathop{ycos}(t)\right) \]
| (%i18) | yy(t); |
\[\]\[\tag{%o18} 3.293349942862704 {{0.9486832980505139}^{t}} \cos{\left( 0.567450475278837 t\mathop{+}1.2622833124845299\right) }\]
| (%i19) | wxplot2d(yy(t),[t,0,50], [style, [lines, 3, 3]]); |
\[\]\[\tag{%t19} \]
\[\]\[\tag{%o19} \]
最初の赤色のグラフと同様のグラフが描けた。最初のy(t)と cosの関数形のyy(t)は同一と考えられる。
| --> | /* 念のためにy(t)とyy(t)の差異を調べてみる */ |
| (%i20) | ydif(t):=y(t)−yy(t); |
\[\]\[\tag{%o20} \mathop{ydif}(t)\mathop{:=}\mathop{y}(t)\mathop{-}\mathop{yy}(t)\]
| (%i21) | wxplot2d(ydif(t),[t,0,50], [style, [lines, 4, 1]]); |
\[\]\[\tag{%t21} \]
\[\]\[\tag{%o21} \]
| --> | /* y(t)とyy(t)の差異は一見すると大きいように見えるが小数点以下15桁の微小な差異である。 */ |
| (%i22) | ydif(t):=round((y(t)−yy(t))·1000)/1000; |
\[\]\[\tag{%o22} \mathop{ydif}(t)\mathop{:=}\frac{\mathop{round}\left( \left( \mathop{y}(t)\mathop{-}\mathop{yy}(t)\right) 1000\right) }{1000}\]
| (%i23) | wxplot2d(ydif(t),[t,0,50], [style, [lines, 4, 1]],[y,−0.0001,0.0001]); |
\[\]\[\tag{%t23} \]
\[\]\[\tag{%o23} \]
| --> | /* y(t)とyy(t)の差異は実務的にはゼロと見なせるだろう。 */ |
\[\]\[\tag{%o20} \]
参考文献
ウォルター・エンダース著 新谷元嗣、藪友良 訳(2019)実証のための計量時系列分析
有斐閣(Walter Enders Applied Econometric time series
Wiley Sons Inc)
Walter Enders 1995 Applied Econometric time series
Wiley Sons Inc
小山昭雄 1995 経済数学教室７ ダイナミック・システム（上） 岩波書店
Enders Applied econometric time series の計算例を参考にしている
;
ARMA(1,1)モデルの反転可能性、因果性等をMAXIMAで探究 invertibility,causality of ARMA(1,1) model with MAXIMAE
Enders Applied econometric time series の計算例を参考にしている