[[FrontPage]]
#contents
2011/06/20からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/10/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* シミュレーション [#k053556f]
シミュレーションは、自然現象や社会現象をモデルを使って再現するもので、
コンピュータ上で専用のアプリケーションを使って行われています。
ここでは、Sageを使って参考文献１「MathematicaによるOR」4章で紹介されている経済モデル
をシミュレーションする方法を紹介します。
モデルの作成には、論理的なモデルを使用する場合とデータに適合する数学モデルを
データから見つける方法があります。数学モデルの場合には、モデルが実際データを
うまく表現しているかを調べるためにデータを教師用データと検証用データに分けて
モデルの表現能力を確認します。
** 経済データ [#b0efec6d]
モデルを作成する前に、今回使用する経済データがどのようなものか見てみましょう。
経済データには、1965年から1989年までの、所得、消費、投資、政府支出が含まれて
います。
それでは、Sageを使って経済データを読み込み、どのような分布をしているか調べてみましょう。
*** データの読み込み [#d5382cb3]
sageでは、豊富なPythonのライブラリを使うことができるため、
データの読み込み処理がとても楽に行えます。
今回は、csvライブラリを使ってカンマ区切りのCSV形式のデータ
(Econometrics.txt)を読み込み、
Year, Ct, Yt, It, Gtの値をリスト変数にセットします。
sageへの入力:
#pre{{
# テストデータ（CSV形式）の読み込み
import csv
fileName = DATA + 'Econometrics.txt'
csvfile = open(fileName) # CSVファイルの読み込み
reader = csv.reader(csvfile) # CSVリーダーの取得
hdr = reader.next() # ヘッダを読む
# データを読む
Year = []; CtData =[]; YtData = []; ItData = []; GtData = [];
for row in reader:
Year.append(N(row[0]))
CtData.append(N(row[1]))
YtData.append(N(row[2]))
ItData.append(N(row[3]))
GtData.append(N(row[4]))
csvfile.close()
# print hdr
# print Year, YtData, CtData, ItData, GtData
}}
*** データのプロット [#q29ddda3]
以下に読み込んだ政府支出、消費、投資、所得をプロットしています。
これから、消費と所得の類似性や投資の変動が大きいことが見て取れます。
sageへの入力:
#pre{{
# データの関係
gtPlot = list_plot(GtData, rgbcolor='red', figsize=3)
ctPlot = list_plot(CtData, rgbcolor='red', figsize=3)
itPlot = list_plot(ItData, rgbcolor='red', figsize=3)
ytPlot = list_plot(YtData, rgbcolor='red', figsize=3)
# プロット表示
html.table([["政府支出", "消費"], [gtPlot, ctPlot], ["投資", "所得"], [ itPlot, ytPlot]])
}}
&ref(out1.png);
** モデルの作成 [#n8865b7e]
経済データを表現するモデルとして以下のような関係をもつ数学モデルを作ります。
&ref(model.png);
矩形はモデルの変数を表し、
Y:所得、C:消費、I:投資、G:政府支出です。
矢印の方向は影響を与える向きを表します。
消費は自分自身へ矢印があり、これは前年度の消費が今年度の
消費に影響を与えていることを表しています。
このモデルを連立方程式で表すと、以下のようになります。
$$
\left\{ \begin{array}{l c l}
C_t & = & a_0 + a_1 Y_t + a_2 C_{t-1} \\
I_t & = & b_0 + b_1 Y_t \\
Y_t & = & C_t + I_t + G_t \\
\end{array} \right.
$$
モデルの作成手順は、以下のようにします。
+ 所得、投資に含まれる係数$a_0, a_1, a_2, b_0, b_1$をデータから求める
+ 連立方程式を解いて、消費、投資、所得を算出する関数を求める
** 係数の決定 [#y61852f9]
連立方程式のCt, Itの係数$a_0, a_1, a_2, b_0, b_1$の値を回帰分析
(find_fit関数)を使って求めます。
教師用データとして、CtList, YtList, GtListには、1966年から1984年のデータをセットし、
Ct1Listには、1年ずらしたCtの値をセットします。
sageへの入力:
#pre{{
# Ctの計算には前年度のCtであるCt-1を使うため1966年から1984年のデータを作成
CtList = CtData[1:20]
Ct1List= CtData[0:19]
YtList = YtData[1:20]
GtList = GtData[1:20]
ItList = ItData[1:20]
}}
*** 消費の係数を求める [#w25581b5]
最初に、消費Ctの式$C_t = a_0 + a_1 Y_t + a_2 C_{t-1}$
の係数$a_0, a_1, a_2$を求めてみましょう。
find_fit関数の引数は以下のように与えます。
#pre{{
find_fit(data, model, options)
data: 変数1の値, 変数2の値, ... , 変数nの値, 関数の値をタプルとするリストをセット
model: model(変数1の値, 変数2の値, ... , 変数nの値)の引数を持つ関数をセット
options: solution_dict=Trueを指定すると係数名をキーとする辞書が戻されます
}}
最初に、変数Yt, Ct, Ct1と係数a0, a1, a2を定義します。
次に、消費のモデルとしてCtModelシンボリック関数を定義します。
データの形式は、消費のモデルCtModelの引数がYt, Ct1であり、
モデルの表す値がCtですから、Yt, Ct1, Ctのタプルのリストを渡します。
各年度のYt, Ct1, Ctのタプルは、zip関数を使って作成します。
sageへの入力:
#pre{{
# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成
data = zip(YtList, Ct1List, CtList)
# モデル関数と変数を定義
(Yt, Ct, Ct1, a0, a1, a2) = var('Yt Ct Ct1 a0 a1 a2')
CtModel(Yt, Ct1) = a0 + a1*Yt + a2*Ct1
# 最適な解のa0, a1, a2を求める
CtFit = find_fit(data, CtModel, solution_dict=True); CtFit
}}
sageの出力:
#pre{{
{a2: 0.59268191617908039, a1: 0.21857976971146675, a0: 7518.3096037374444}
}}
*** 投資の係数を求める [#q20e3942]
消費と同様に、投資Itの式$I_t = b_0 + b_1 Y_t$
の係数$b_0, b_1$を求めます。
新たに変数Yt, b0, b1を定義し、投資のモデルとしてItModelシンボリック関数を定義します。
データの形式は、投資のモデルItModelの引数がYtであり、
モデルの表す値がItですから、Yt, Itのタプルのリストを渡します。
sageへの入力:
#pre{{
# 同様にItに対するb0, b1を求める
data = zip(YtList, ItList)
(Yt, b0, b1) = var('Yt b0 b1')
ItModel(Yt) = b0 + b1*Yt
ItFit = find_fit(data, ItModel, solution_dict=True); ItFit
}}
sageの出力:
#pre{{
{b1: 0.15818655082721944, b0: 311.21923398160885}
}}
** 連立方程式を解く [#w1a41a89]
係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連立方程式
を解きます。
係数を入力して式を定義する代わりに、
#pre{{
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
}}
のようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1, a2を代入しています。
これで、モデルを変更しても再計算が楽にできます。
解の結果、Ct, It, YtはGtとCt1からなる式に変換されています。
経済モデルから前年度の消費と今年度の政府支出によって、今年度の消費、投資、所得を求める式を得ることができました。
sageへの入力:
#pre{{
# 連立方程式をCt, It, Ytに対して解く
(Gt, It) = var('Gt It')
CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)
ItEq = (It == ItModel(Yt)).subs(ItFit)
YtEq = (Yt == Ct + It + Gt)
eqn = [CtEq, ItEq, YtEq]; eqn
}}
sageの出力:
#pre{{
[Ct == 0.59268191617908039*Ct1 + 0.21857976971146675*Yt + 7518.3096037374444,
It == 0.15818655082721944*Yt + 311.21923398160885, Yt == Ct + Gt + It]
}}
sageへの入力:
#pre{{
sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol
}}
sageの出力:
#pre{{
[{Ct: 327042806728/408524384645*Ct1 + 59574711/169864609*Gt +
9519709947510127/927460765140, Yt: 388497958456/408524384645*Ct1 +
272553640/169864609*Gt + 582572871301814/46373038257,
It: 61455151728/408524384645*Ct1 + 43114320/169864609*Gt +
710582492842051/309153588380}]
}}
** 計算結果と実データの比較 [#n723b432]
計算結果と実データを比較するために、zip関数を使って前年度消費と政府支出をタプルとするリストを作成します。
sageへの入力:
#pre{{
# 入力データをセット
data = zip(Ct1List, GtList)
}}
*** 消費・投資・所得計算関数の定義 [#a82eb912]
連立方程式から消費Ct、投資It、所得Ytを計算する式を取り出し、
前年度消費Ct1と今年度政府支出を引数とする関数ctFunc, itFunc, ytFuncを定義します。
sageへの入力:
#pre{{
ctFunc(Ct1, Gt) = sol[0][Ct]; view(ctFunc) # 消費計算関数
itFunc(Ct1, Gt) = sol[0][It]; view(itFunc) # 投資計算関数
ytFunc(Ct1, Gt) = sol[0][Yt]; view(ytFunc) # 所得計算関数
}}
&ref(out2.png);
*** 教師用データに対するシミュレーション [#q2829564]
求めた関数ctFunc, itFunc, ytFuncを使って教師用データの消費、投資、所得をシミュレーション
してみましょう。
dataの各タプルに対し、リスト内包を使って関数から計算結果のリストを作成します。
sageへの入力:
#pre{{
# 消費のシミュレーション結果
calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]
# 投資のシミュレーション結果
calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data]
# 所得のシミュレーション結果
calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]
}}
*** 教師用データに対するシミュレーション結果 [#ya7a61c8]
list_plotを使って教師用データに対するシミュレーション結果（青い線）と実データ（赤い点）をプロットします。
単純なモデルの割には消費と所得に対しては、シミュレーション結果と実データが合っています。
投資については、所得Ytのみが影響するモデルとなっているため、実データの変動をうまく表現できていません。
sageへの入力:
#pre{{
# 計算値（青）と実測値（赤）のプロット
gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)
calCtPlot = list_plot(calCtList, plotjoined=true, figsize=3)
ctPlot = list_plot(CtList, rgbcolor='red', figsize=3)
calItPlot = list_plot(calItList, plotjoined=true, figsize=3)
itPlot = list_plot(ItList, rgbcolor='red', figsize=3)
calYtPlot = list_plot(calYtList, plotjoined=true, figsize=3)
ytPlot = list_plot(YtList, rgbcolor='red', figsize=3)
# プロット表示
html.table([["政府支出", "消費"], [gtPlot, (calCtPlot + ctPlot)], ["投資", "所得"], [(calItPlot + itPlot), (calYtPlot + ytPlot)]])
}}
&ref(out3.png);
** 検証用データとの比較 [#lfa4a75b]
同様に検証用データに対してシミュレーションを行ってみましょう。
消費、所得については、モデルはよく表現できていますが、
投資については急激な増加をうまく表現できていません。
sageへの入力:
#pre{{
# 検証用データ
valCtList = [calCtList[18]]
valItList = [calItList[18]]
valYtList = [calYtList[18]]
Ct1List = CtData[0:24]
GtList = GtData[1:25]
CtList = CtData[1:25]
for i in range(19, 24):
valCtList.append(ctFunc(Ct1List[i], GtList[i]))
valItList.append(itFunc(Ct1List[i], GtList[i]))
valYtList.append(ytFunc(Ct1List[i], GtList[i]))
}}
sageへの入力:
#pre{{
# 教師用＋検証用の実測値（赤）
ctPlot = list_plot(CtData[1:24], rgbcolor='red', figsize=3)
itPlot = list_plot(ItData[1:25], rgbcolor='red', figsize=3)
ytPlot = list_plot(YtData[1:25], rgbcolor='red', figsize=3)
gtPlot = list_plot(GtData[1:25], rgbcolor='red', figsize=3)
# 計算値（青）、検証値（灰）
valCtPlot = list_plot(zip(range(18,24), valCtList), plotjoined=true, rgbcolor ='gray', figsize=3)
valItPlot = list_plot(zip(range(18,24), valItList), plotjoined=true, rgbcolor ='gray', figsize=3)
valYtPlot = list_plot(zip(range(18,24), valYtList), plotjoined=true, rgbcolor ='gray', figsize=3)
# プロット
ctPlt = (calCtPlot + ctPlot +valCtPlot)
itPlt = (calItPlot + itPlot + valItPlot)
ytPlt = (calYtPlot + ytPlot + valYtPlot)
html.table([["政府支出", "消費"], [gtPlot, ctPlt], ["投資", "所得"], [itPlt, ytPlt]])
}}
&ref(out4.png);
** 予測 [#u69fd97a]
シミュレーションの結果から、モデルは消費をうまく表現していることが確認でます。
このことから、もし政府支出が予測できたら今年度度消費の値から来年度の消費、投資、所得
を計算することができます。
問題を簡単にするために、政府支出Gtを$G_t(x) = c_0 + c_1 x$として回帰分析をします。
sageへの入力:
#pre{{
# 1985年から1989年のデータを予測
var('x c0 c1')
data = zip(range(20), GtList[0:20])
# 直線回帰から政府支出のモデルGtの係数を求める
GtModel(x) = c0 + c1*x
GtFit = find_fit(data, GtModel, solution_dict=True); GtFit
}}
sageの出力:
#pre{{
{c1: 1837.6999989361932, c0: 36476.850068506872}
}}
*** 政府支出の予測結果 [#w03f99d7]
直線回帰で求めた政府支出と実データをプロットすると以下のようになります。
sageへの入力:
#pre{{
gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)
gtFunc(x) = GtModel.subs(GtFit)
calGtPlot = plot(gtFunc, [x, 0, 19], figsize=3)
(calGtPlot + gtPlot).show()
}}
&ref(out5.png);
*** 1967年から1989年までの予測 [#g757862e]
ちょっと乱暴な遊びですが、
政府支出関数gtFuncを使って、1967年から1989年まで政府支出を求め、前年度の消費Ct1List
の値とgtFuncの値から消費、投資、所得を順番に計算してみましょう。
sageへの入力:
#pre{{
# 予測
predCt1List =[ Ct1List[0]] + range(23)
predCtList = [calCtList[0]]
predItList = [calItList[0]]
predYtList = [calYtList[0]]
for i in range(1, 24):
predCt1List[i] = ctFunc(predCt1List[i-1], gtFunc(i-1))
predCtList.append(ctFunc(predCt1List[i], gtFunc(i)))
predItList.append(itFunc(predCt1List[i], gtFunc(i)))
predYtList.append(ytFunc(predCt1List[i], gtFunc(i)))
}}
*** 予測結果 [#a7249c61]
政府支出の乱暴な予測にも関わらず、予測値(緑)は先のシミュレーション結果と
ほぼ同じ傾向を示しています。
このようにシミュレーションは、未来の予測にも活用できることが理解して頂けたでしょうか。
sageへの入力:
#pre{{
# 計算値（青）、予測値（緑）と実測値（赤）
predCtPlot = list_plot(zip(range(24), predCtList), plotjoined=true, rgbcolor ='green', linestyle='dashed', figsize=3)
predItPlot = list_plot(zip(range(24), predItList), plotjoined=true, rgbcolor ='green', linestyle='dashed', figsize=3)
predYtPlot = list_plot(zip(range(24), predYtList), plotjoined=true, rgbcolor ='green', linestyle='dashed', figsize=3)
calGtPlot = plot(gtFunc, [x, 0, 23], figsize=3)
# プロット
ctPlt = (calCtPlot + predCtPlot + ctPlot + valCtPlot)
itPlt = (calItPlot + predItPlot + itPlot + valItPlot)
ytPlt = (calYtPlot + predYtPlot + ytPlot + valYtPlot)
html.table([["政府支出", "消費"], [(calGtPlot + gtPlot), ctPlt], ["投資", "所得"], [itPlt, ytPlt]])
}}
&ref(out6.png);
** コメント [#e7a9ce66]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
#vote(おもしろかった[3],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
