{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "\n", "\tHiroshi TAKEMOTO\n", "\t(take.pwave@gmail.com)\n", "\t

シミュレーション

\n", "\t

\n", "\t\tシミュレーションは、自然現象や社会現象をモデルを使って再現するもので、\n", "\t\tコンピュータ上で専用のアプリケーションを使って行われています。\n", "\t\tここでは、Sageを使って「MathematicaによるOR」4章で紹介されている経済モデル\n", "\t\tをシミュレーションする方法を紹介します。\n", "\t

\n", "\t

\n", "\t\tモデルの作成には、論理的なモデルを使用する場合とデータに適合する数学モデルを\n", "\t\tデータから見つける方法があります。数学モデルの場合には、モデルが実際データを\n", "\t\tうまく表現しているかを調べるためにデータを教師用データと検証用データに分けて\n", "\t\tモデルの表現能力を確認します。\t\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%HTML\n", "" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# jupyter用のdisplayメソッド\n", "from IPython.display import display, Latex, HTML, Math, JSON\n", "# sageユーティリティ\n", "load('script/sage_util.py')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

経済データ

\n", "\t

\n", "\t\tモデルを作成する前に、今回使用する経済データがどのようなものか見てみましょう。\n", "\t\t経済データには、1965年から1989年までの、所得、消費、投資、政府支出が含まれて\n", "\t\tいます。\n", "\t

\n", "\t

\n", "\t\tそれでは、Sageを使って経済データを読み込み、どのような分布をしているか調べてみましょう。\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

データの読み込み

\n", "\t

\n", "\t\tsageでは、豊富なPythonのライブラリを使うことができるため、\n", "\t\tデータの読み込み処理がとても楽に行えます。\n", "\t

\n", "\t

\n", "\t\t今回は、csvライブラリを使ってカンマ区切りのCSV形式のデータ\n", "\t\t(Econometrics.txt)を読み込み、\n", "\t\tYear, Ct, Yt, It, Gtの値をリスト変数にセットします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# テストデータ(CSV形式)の読み込み\n", "import csv\n", "fileName = 'data/Econometrics.txt'\n", "csvfile = open(fileName) # CSVファイルの読み込み\n", "reader = csv.reader(csvfile) # CSVリーダーの取得\n", "hdr = reader.next() # ヘッダを読む\n", "# データを読む\n", "Year = []; CtData =[]; YtData = []; ItData = []; GtData = []; \n", "for row in reader:\n", " Year.append(N(row[0]))\n", " CtData.append(N(row[1]))\n", " YtData.append(N(row[2]))\n", " ItData.append(N(row[3]))\n", " GtData.append(N(row[4]))\n", "csvfile.close()\n", "# print hdr\n", "# print Year, YtData, CtData, ItData, GtData" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

データのプロット

\n", "\t

\n", "\t\t以下に読み込んだ政府支出、消費、投資、所得をプロットしています。\n", "\t\tこれから、消費と所得の類似性や投資の変動が大きいことが見て取れます。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
政府支出消費
投資所得
" ], "text/plain": [ "<__main__.Table2Html instance at 0x7fa32d978b48>" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# データの関係\n", "gtPlot = list_plot(GtData, rgbcolor='red', figsize=3)\n", "ctPlot = list_plot(CtData, rgbcolor='red', figsize=3)\n", "itPlot = list_plot(ItData, rgbcolor='red', figsize=3)\n", "ytPlot = list_plot(YtData, rgbcolor='red', figsize=3)\n", "# プロット表示\n", "Table2Html([[\"政府支出\", \"消費\"], [_to_png(gtPlot), _to_png(ctPlot)], \n", " [\"投資\", \"所得\"], [ _to_png(itPlot), _to_png(ytPlot)]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

モデルの作成

\n", "\t

\n", "\t\t経済データを表現するモデルとして以下のような関係をもつ数学モデルを作ります。\n", "\t

\n", "\t

\n", "\t

\n", "\t\t矩形はモデルの変数を表し、\n", "\t\tY:所得、C:消費、I:投資、G:政府支出です。\n", "\t\t矢印の方向は影響を与える向きを表します。\n", "\t\t消費は自分自身へ矢印があり、これは前年度の消費が今年度の\n", "\t\t消費に影響を与えていることを表しています。\t\t\n", "\t

\n", "\t

\n", "\t\tこのモデルを連立方程式で表すと、以下のようになります。\n", "$$\n", "\t\\left\\{ \\begin{array}{l c l}\n", "\t\tC_t & = & a_0 + a_1 Y_t + a_2 C_{t-1} \\\\\n", "\t\tI_t & = & b_0 + b_1 Y_t \\\\\n", "\t\tY_t & = & C_t + I_t + G_t \\\\\n", "\t\\end{array} \\right.\n", "$$\t\t\n", "\t

\n", "\t

\n", "\t\tモデルの作成手順は、以下のようにします。\n", "\t\t

    \n", "\t\t\t
  1. 所得、投資に含まれる係数$a_0, a_1, a_2, b_0, b_1$をデータから求める
  2. \n", "\t\t\t
  3. 連立方程式を解いて、消費、投資、所得を算出する関数を求める
  4. \n", "\t\t
\n", "\t

\n", "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

係数の決定

\n", "\t

\n", "\t\t連立方程式のCt, Itの係数$a_0, a_1, a_2, b_0, b_1$の値を回帰分析\n", "\t\t(find_fit関数)を使って求めます。\t\t\n", "\t

\n", "\t

\n", "\t\t教師用データとして、CtList, YtList, GtListには、1966年から1984年のデータをセットし、\n", "\t\tCt1Listには、1年ずらしたCtの値をセットします。\t\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Ctの計算には前年度のCtであるCt-1を使うため1966年から1984年のデータを作成\n", "CtList = CtData[1:20]\n", "Ct1List= CtData[0:19]\n", "YtList = YtData[1:20]\n", "GtList = GtData[1:20]\n", "ItList = ItData[1:20]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

消費の係数を求める

\n", "\t

\n", "\t\t最初に、消費Ctの式$C_t = a_0 + a_1 Y_t + a_2 C_{t-1}$\n", "\t\tの係数$a_0, a_1, a_2$を求めてみましょう。\n", "\t

\n", "\t

\n", "\t\tfind_fit関数の引数は以下のように与えます。\n", "\t\t

\n",
    "find_fit(data, model, options)\n",
    "data: 変数1の値, 変数2の値, ... , 変数nの値, 関数の値をタプルとするリストをセット\n",
    "model: model(変数1の値, 変数2の値, ... , 変数nの値)の引数を持つ関数をセット\n",
    "options: solution_dict=Trueを指定すると係数名をキーとする辞書が戻されます\n",
    "\t\t
\n", "\t

\n", "\t

\n", "\t\t最初に、変数Yt, Ct, Ct1と係数a0, a1, a2を定義します。\n", "\t\t次に、消費のモデルとしてCtModelシンボリック関数を定義します。\n", "\t

\n", "\t

\n", "\t\tデータの形式は、消費のモデルCtModelの引数がYt, Ct1であり、\n", "\t\tモデルの表す値がCtですから、Yt, Ct1, Ctのタプルのリストを渡します。\n", "\t

\n", "\t

\n", "\t\t各年度のYt, Ct1, Ctのタプルは、zip関数を使って作成します。\n", "\t

\n", "\t\t\n", "" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{a1: 0.21857976971146675, a2: 0.5926819161790804, a0: 7518.309603737444}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 各年度のYt, Ct1, Ctをタプルにまとめたリストを作成\n", "data = zip(YtList, Ct1List, CtList)\n", "# モデル関数と変数を定義\n", "(Yt, Ct, Ct1, a0, a1, a2) = var('Yt Ct Ct1 a0 a1 a2')\n", "CtModel(Yt, Ct1) = a0 + a1*Yt + a2*Ct1\n", "# 最適な解のa0, a1, a2を求める\n", "CtFit = find_fit(data, CtModel, solution_dict=True); CtFit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

投資の係数を求める

\n", "\t

\n", "\t\t消費と同様に、投資Itの式$I_t = b_0 + b_1 Y_t$\n", "\t\tの係数$b_0, b_1$を求めます。\n", "\t

\n", "\t

\n", "\t\t新たに変数Yt, b0, b1を定義し、投資のモデルとしてItModelシンボリック関数を定義します。\n", "\t

\n", "\t

\n", "\t\tデータの形式は、投資のモデルItModelの引数がYtであり、\n", "\t\tモデルの表す値がItですから、Yt, Itのタプルのリストを渡します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{b0: 311.21923398160885, b1: 0.15818655082721944}" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 同様にItに対するb0, b1を求める\n", "data = zip(YtList, ItList)\n", "(Yt, b0, b1) = var('Yt b0 b1')\n", "ItModel(Yt) = b0 + b1*Yt\n", "ItFit = find_fit(data, ItModel, solution_dict=True); ItFit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

連立方程式を解く

\n", "\t

\n", "\t\t係数a0, a1, b2, b0, b1が決まったので、solve関数を使って連立方程式\n", "\t\tを解きます。\n", "\t

\n", "\t

\n", "\t\t係数を入力して式を定義する代わりに、\n", "

\n",
    "CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)\n",
    "
\n", "\t\tのようにsubs関数を使ってCt == a0 + a1*Yt + a2*Ct1のa0, a1, a2を代入しています。\n", "\t\tこれで、モデルを変更しても再計算が楽にできます。\t\n", "\t\t解の結果、Ct, It, YtはGtとCt1からなる式に変換されています。\t\n", "\t

\n", "\t

\n", "\t\t経済モデルから前年度の消費と今年度の政府支出によって、今年度の消費、投資、所得を求める式を得ることができました。\n", "\t

" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[Ct == 0.5926819161790804*Ct1 + 0.21857976971146675*Yt + 7518.309603737444,\n", " It == 0.15818655082721944*Yt + 311.21923398160885,\n", " Yt == Ct + Gt + It]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 連立方程式をCt, It, Ytに対して解く\n", "(Gt, It) = var('Gt It')\n", "CtEq = (Ct == CtModel(Yt, Ct1)).subs(CtFit)\n", "ItEq = (It == ItModel(Yt)).subs(ItFit)\n", "YtEq = (Yt == Ct + It + Gt)\n", "eqn = [CtEq, ItEq, YtEq]; eqn" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "[{Yt: 1760421695327231515898/1851168494654738992935*Ct1 + 1665464251735955/1037973413620683*Gt + 27161739944082251799229948075/2162085545220791309256249,\n", " Ct: 7409733296710051049024/9255842473273694964675*Ct1 + 364036792607125/1037973413620683*Gt + 22192235380675639267030739855/2162085545220791309256249,\n", " It: 1392375179926106530466/9255842473273694964675*Ct1 + 263454045508147/1037973413620683*Gt + 4969504563406612532199208220/2162085545220791309256249}]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "sol = solve(eqn, [Ct, It, Yt], solution_dict=True); sol" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

計算結果と実データの比較

\n", "\t

\n", "\t\t計算結果と実データを比較するために、zip関数を使って前年度消費と政府支出をタプルとするリストを作成します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 入力データをセット\n", "data = zip(Ct1List, GtList)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

消費・投資・所得計算関数の定義

\n", "\t

\n", "\t\t連立方程式から消費Ct、投資It、所得Ytを計算する式を取り出し、\n", "\t\t前年度消費Ct1と今年度政府支出を引数とする関数ctFunc, itFunc, ytFuncを定義します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "(Ct1, Gt) |--> 7409733296710051049024/9255842473273694964675*Ct1 + 364036792607125/1037973413620683*Gt + 22192235380675639267030739855/2162085545220791309256249" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "(Ct1, Gt) |--> 1392375179926106530466/9255842473273694964675*Ct1 + 263454045508147/1037973413620683*Gt + 4969504563406612532199208220/2162085545220791309256249" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "(Ct1, Gt) |--> 1760421695327231515898/1851168494654738992935*Ct1 + 1665464251735955/1037973413620683*Gt + 27161739944082251799229948075/2162085545220791309256249" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ctFunc(Ct1, Gt) = sol[0][Ct]; show(ctFunc ) # 消費計算関数\n", "itFunc(Ct1, Gt) = sol[0][It]; show(itFunc) # 投資計算関数\n", "ytFunc(Ct1, Gt) = sol[0][Yt]; show(ytFunc) # 所得計算関数" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

教師用データに対するシミュレーション

\n", "\t

\n", "\t\t求めた関数ctFunc, itFunc, ytFuncを使って教師用データの消費、投資、所得をシミュレーション\n", "\t\tしてみましょう。\n", "\t

\n", "\t

\n", "\t\tdataの各タプルに対し、リスト内包を使って関数から計算結果のリストを作成します。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 消費のシミュレーション結果\n", "calCtList = [ctFunc(Ct1, Gt).n() for (Ct1, Gt) in data]\n", "# 投資のシミュレーション結果\n", "calItList = [itFunc(Ct1, Gt) for (Ct1, Gt) in data]\n", "# 所得のシミュレーション結果\n", "calYtList = [ytFunc(Ct1, Gt) for (Ct1, Gt) in data]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

教師用データに対するシミュレーション結果

\n", "\t

\n", "\t\tlist_plotを使って教師用データに対するシミュレーション結果(青い線)と実データ(赤い点)をプロットします。\n", "\t\t単純なモデルの割には消費と所得に対しては、シミュレーション結果と実データが合っています。\n", "\t\t投資については、所得Ytのみが影響するモデルとなっているため、実データの変動をうまく表現できていません。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
政府支出消費
投資所得
" ], "text/plain": [ "<__main__.Table2Html instance at 0x7fa32da20128>" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 計算値(青)と実測値(赤)のプロット\n", "gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)\n", "calCtPlot = list_plot(calCtList, plotjoined=true, figsize=3)\n", "ctPlot = list_plot(CtList, rgbcolor='red', figsize=3)\n", "calItPlot = list_plot(calItList, plotjoined=true, figsize=3)\n", "itPlot = list_plot(ItList, rgbcolor='red', figsize=3)\n", "calYtPlot = list_plot(calYtList, plotjoined=true, figsize=3)\n", "ytPlot = list_plot(YtList, rgbcolor='red', figsize=3)\n", "# プロット表示\n", "Table2Html([[\"政府支出\", \"消費\"], [_to_png(gtPlot), _to_png(calCtPlot + ctPlot)], \n", " [\"投資\", \"所得\"], [_to_png(calItPlot + itPlot), _to_png(calYtPlot + ytPlot)]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

検証用データとの比較

\n", "\t

\n", "\t\t同様に検証用データに対してシミュレーションを行ってみましょう。\n", "\t

\n", "\t

\n", "\t\t消費、所得については、モデルはよく表現できていますが、\n", "\t\t投資については急激な増加をうまく表現できていません。\t\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 検証用データ\n", "valCtList = [calCtList[18]]\n", "valItList = [calItList[18]]\n", "valYtList = [calYtList[18]]\n", "Ct1List = CtData[0:24]\n", "GtList = GtData[1:25]\n", "CtList = CtData[1:25]\n", "for i in range(19, 24):\n", " valCtList.append(ctFunc(Ct1List[i], GtList[i]))\n", " valItList.append(itFunc(Ct1List[i], GtList[i]))\n", " valYtList.append(ytFunc(Ct1List[i], GtList[i]))" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "
政府支出消費
投資所得
" ], "text/plain": [ "<__main__.Table2Html instance at 0x7fa3263f3290>" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 教師用+検証用の実測値(赤)\n", "ctPlot = list_plot(CtData[1:24], rgbcolor='red', figsize=3)\n", "itPlot = list_plot(ItData[1:25], rgbcolor='red', figsize=3)\n", "ytPlot = list_plot(YtData[1:25], rgbcolor='red', figsize=3)\n", "gtPlot = list_plot(GtData[1:25], rgbcolor='red', figsize=3)\n", "# 計算値(青)、検証値(灰)\n", "valCtPlot = list_plot(zip(range(18,24), valCtList), plotjoined=true, rgbcolor ='gray', figsize=3)\n", "valItPlot = list_plot(zip(range(18,24), valItList), plotjoined=true, rgbcolor ='gray', figsize=3)\n", "valYtPlot = list_plot(zip(range(18,24), valYtList), plotjoined=true, rgbcolor ='gray', figsize=3)\n", "# プロット\n", "ctPlt = (calCtPlot + ctPlot +valCtPlot)\n", "itPlt = (calItPlot + itPlot + valItPlot)\n", "ytPlt = (calYtPlot + ytPlot + valYtPlot)\n", "Table2Html([[\"政府支出\", \"消費\"], [_to_png(gtPlot), _to_png(ctPlt)], \n", " [\"投資\", \"所得\"], [_to_png(itPlt), _to_png(ytPlt)]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

予測

\n", "\t

\n", "\t\tシミュレーションの結果から、モデルは消費をうまく表現していることが確認でます。\n", "\t\tこのことから、もし政府支出が予測できたら今年度度消費の値から来年度の消費、投資、所得\n", "\t\tを計算することができます。\n", "\t

\n", "\t

\n", "\t\t問題を簡単にするために、政府支出Gtを$G_t(x) = c_0 + c_1 x$として回帰分析をします。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "{c0: 36476.85006850687, c1: 1837.6999989361932}" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 1985年から1989年のデータを予測\n", "var('x c0 c1')\n", "data = zip(range(20), GtList[0:20])\n", "# 直線回帰から政府支出のモデルGtの係数を求める\n", "GtModel(x) = c0 + c1*x\n", "GtFit = find_fit(data, GtModel, solution_dict=True); GtFit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

政府支出の予測結果

\n", "\t

\n", "\t\t直線回帰で求めた政府支出と実データをプロットすると以下のようになります。\n", "\t

\t\n", "" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAARwAAADSCAYAAACYYX+QAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3XtcVHX6B/DPACGCchlxZrwk4gW0QCI2L+S6JshFNllX\nakXsqluCWParcFNTy9qk3Wovathmq/UzwW2zsBRNAlGWa15AwdLEH7Y1RwkCJeUyPL8/vnL0OKAM\nMjPAPO/Xa17NmfNlzjMTPnzP96oiIgJjjFmAnbUDYIzZDk44jDGL4YTDGLMYTjiMMYvhhMMYsxhO\nOIwxi+GEwxizGE44jDGLsamEQ0Soq6sDj3VkzDpsKuFcuHABbm5uuHDhgrVDYaxnOnYMePxxYNEi\n4IcfTP5xBzOExBjrjaqqgF/9CqiuFsf79wOlpYBK1eG3sKkaTqs5c+Zg5syZ2LZtm7VDYaznOH78\narJp67gDVLY0ebOurg5ubm6ora2Fq6urtcNhrGeRJMDXF6itFcejRwMnTgB2Ha+32GQNhzHWhrfe\nEklk8mSRSK6n1QL79gG/+x3w6KPiuQnJBuAaDmMMAHJyRPtMq4AA4MiRLr8M13AYY0BFxY2Puwgn\nHMYYMH06oNFcPZ471yyXsclbqsjISDg4OCA2NhaxsbHWDoux7qGyEvjoI0CnA2JjTeru7iibTDjc\nhsOYdfAtFWO2oqYGMBisGgInHMZ6u59/Bu67D1CrgaFDgcOHrRYKJxzGeruNG4HsbPFcrweeecZq\noXDCYay3+/ln5XF9vXXigI0mHJ5LxWzK448Dw4eL546OwIsvWi0U7qVizBbU1QGHDonE05p8rICX\np2DMFri6AlOnWjsK27ylYoxZByccxpjFcMJhjFmMTSYc7qVizDpMSjje3t6ws7MzeixevFgus3Ll\nSgwePBjOzs6YPn06Tp06pXiPhoYGLFq0CJ6enujfvz9iYmJw7tw5RZmamhrExcXBzc0NHh4eWLBg\nAeqvGztw9uxZREVFwcXFBTqdDklJSWhpaenQ50hNTUV6ejpP3GTd3vnzQHIy0Gv6kskEVVVVJEmS\n/Ni3bx/Z2dlRTk4OERGtXbuWPDw8aOfOnVRaWkrR0dE0YsQIamhokN9j4cKF5OXlRdnZ2XTo0CGa\nNGkSTZ48WXGdiIgICgwMpKKiIsrNzaXRo0dTXFycfN5gMJCfnx+FhYVRSUkJZWRk0MCBA2n58uU3\njL+2tpYAUG1trSkfmzGLa2kh2rSJSK0Wj1OnrB1R1zAp4Vzv6aefptGjR8vHgwYNojfffFM+rq2t\nJScnJ0pLS5OPHR0d6eOPP5bLnDhxglQqFRUUFBARUVlZGalUKjp06JBcJiMjg+zt7emHH34gIqJd\nu3aRg4MDnT9/Xi6TkpJC7u7u1NTU1G68nHBYT1BWRjRlChFA9PDDROfOWTuirtPpNpympiZs3boV\n8+fPBwBUVFRAr9cjJCRELuPq6ooJEyYgLy8PAFBcXIzm5mZFGV9fXwwbNkwuk5+fDw8PDwQGBspl\nQkNDoVKpUFBQIJfx9/eHp6enXCY8PBy1tbU4fvx4Zz8SY1Z16ZIYBBwQIKY8ZWYCW7YAAwdaO7Ku\n0+mEs2PHDtTW1uKRRx4BAOj1eqhUKmi1WkU5rVYLvV4PAJAkCY6OjkajfK8to9frobl25TEA9vb2\nUKvVijJtXaf1HGM9zb59wLhxwOuvA8uWAUePAtOmWTuqrtfphPPee+8hMjISOp2uK+NhzKacOwfM\nmydW+Bw6FCgpAVavBpycrB2ZeXRqakNlZSX27duHTz75RH5Np9OBiCBJkqL2IUmSfHuk0+nQ2NiI\nuro6RS1HkiQ5cel0OqNeK4PBgOrqakWZoqIiRRlJkuRzNzNnzhw4OCg/Oi83yiyppQXYtAlISgLs\n7cWt00MPmWVVz+6lMw0/q1atosGDB5PBYFC83l6j8fbt2+XjmzUal5eXk52dnaLReM+ePYpG4927\ndxs1Gm/cuJHc3d2psbGx3bi50Zh1B6WlRPfeKxqFH3uM6Jpf4845d44oK4voyr+P7szkhNPS0kJe\nXl60bNkyo3PJycmkVqspPT2dSkpKKDo6mkaNGqXoFo+Pj6fhw4dTVlYWFRcXU3BwsFG3eGRkJAUF\nBVFhYSEdPHiQfHx8aN68efJ5g8FA48aNo4iICDp69ChlZGSQRqOhFStW3DB2TjjMmurriV54gcjB\ngWjMGKLs7C5409JS0W8OELm6EhUWdsGbmo/JCWfv3r1kZ2dHJ0+ebPP8qlWraNCgQdS3b18KCwsz\nKnf58mVKTEykAQMGUL9+/SgmJoYkSVKUqampobi4OHJ1dSV3d3dasGAB1dfXK8pUVlZSVFQUubi4\nkEajoaSkJKMa1/U44TBr2b2byNubqE8fopdfJrp8uYveeP58kWxaHw880EVvbB68Hg5jZtS6omdq\nKhASArz9tthNt8ssXgysW3f1+JFHgM2bu/ACXcsm51IxZm4tLUBKCjBmjBhP88EHwBdfdHGyAYDl\nywE/P/Hcxwd4+eUuvkDXsskaDm+Ex8yptBR44gkgPx9YsEDMhVKrzXhBIrEFjIdHt+/mssmEw7dU\nzBzq60UF4403AF9fsVnC5MnWjqp74SVGGesCu3YBCQmAJImk89xzYr1ypsQJh7Fb8P33wJIlwL/+\nBYSFifaakSOtHVX3xY3GjHWCwQCsXw+MHQvk5AAffghkZHCyuRlOOKxna2kR8wPuuAN44AGgutrs\nlzxyBJg0CUhMBGJjgfJy8d9u3l7bLdjkLVXrXCrupeoF3n0X+NOfxPPycjHr8YMPzHKpixfFxMq/\n/EXUbHJzgeBgs1yq17LJhJOamsq9VL3FdUvY4ttvzXKZnTuBRYuAqirg1VeB//kf4LbbbvFNf/xR\nvOGoUWIGpw3gWyrWfTU1idulKVPEIjEGg3GZ3/xG+S8/Jqbz12tpAYqKgBMn5Je++w6YPRuYORO4\n807g+HFg6dIuSDY7d4r1KMaMAX71K7H6VlvOngWyssQ4m97AmvMqLI3nUvUwq1cr5wm99lrb5QoL\nidasIbpmFQKTNTcTRUbK12pe80f661+J+vUj0umI0tLEOsP00UdEt99ONHQo0ZVVEDpl9GjlZ3v3\nXeMye/cS9e0rzg8eTFRR0fnrdROccFj3NXu28h/l3Lmdf6+jR4mWLSPasEEkl+tlZsrXKcbdFIQi\nUqlaKCGB6KefrpSpqhKzL1vjcXTs/ILDo0bdPOHcd5+yTFJS567VjfAtFeu+IiNvfNxR33wjWnf/\n+EcxOi8hwbhMnz64gH54Bm9iPArRpHJE3sEWrF8PuLldKVNdDTQ0XP2ZxkbRDtMZb7xxdVm/e+8F\n5s41LuPsrDx2cenctboTa2c8S2qt4URGRtL9999PH374obVDYjeTmkr09NPiVqazNmxQ1hS0WqMi\nO3YQDXX5kZxxkf5kl0SNm943fh+DgWjatKvvM3Vq27Wlb78lmjVLlP3ss/bjqqoiKi8nam+nkbIy\ncfsGEE2ceE1Vq+fiuVSs98vKUq5IPnWqeA2iTXbxYuDTT4GoKGD9mmp4jXYE+vVr+70aGoCPPhIp\n54EHgD59jMuMGQN8/bV47ugIHDvW+WniLS1AXR3g7t65n+9mOOEw25CSAvzzn8CQIcC6dWjWDMbf\n/y62ZXFzA/72N+C3v+2CwXuNjcZJaOdO4Ne/vsU37h044TCbU1QEPPmkGDGcmAi88grQpb8O06bJ\nNSgMGCBqOLy7CQAbHfjHbFNdnVivav164K67gIIC4J57zHCh9HTgrbeA2lrg97/nZHMNruGwXo8I\n+Phj4KmnRA5Ys0a02zjwn1uLs8lu8Tlz5mDmzJnYtm2btUNhZnbmDHD//WIA8j33AGVlYo1hTjbW\nYXLC+f777/HQQw/B09MTzs7OCAgIwKFDhxRlVq5cicGDB8PZ2RnTp0/HqevmuzQ0NGDRokXw9PRE\n//79ERMTY7T5XU1NDeLi4uDm5gYPDw8sWLAA9fX1ijJnz55FVFQUXFxcoNPpkJSUhJaWlpt+htTU\nVKSnp/PEzV6sqQn485/FdISjR4EdO4BPPgGGDbN2ZDbOlD70mpoaGj58OM2fP5+Ki4vpzJkz9MUX\nX9Dp06flMmvXriUPDw/auXMnlZaWUnR0NI0YMUKxN9XChQvJy8uLsrOz6dChQzRp0iSjvakiIiIo\nMDCQioqKKDc3l0aPHk1xcXHyeYPBQH5+fhQWFkYlJSWUkZFBAwcOpOXLl7cbP480tg15eUTjxhHZ\n2REtWUJUV2ftiFgrkxLO0qVLacqUKTcs097um2lpafLxzXbfLCsrI5VKpdh9MyMjQ7H75q5du4x2\n30xJSSF3d3dqamcgFSec3u2nn4gSEohUKqKgIKKvvrJ2ROx6Jt1S7dy5E7/4xS/w4IMPQqvV4u67\n78a7774rn6+oqIBer0dISIj8mqurKyZMmIC8vDwAQHFxMZqbmxVlfH19MWzYMLlMfn4+PDw85D3J\nASA0NBQqlQoFBQVyGX9/f3h6esplwsPDUVtbi+PHj5vysZg1/PijGDjn7y8Gw9wCImD7djHe7v33\nxXo1BQXA3Xd3Uaysy5iUcE6fPo23334bvr6+2Lt3L+Lj4/HUU0/hgysLHun1eqhUKmi1WsXPabVa\n6PV6AIAkSXB0dDTqJbq2jF6vh0ajUZy3t7eHWq1WlGnrOq3nWDe3aJEYsXvsmBgI8/77nXqbigpg\nxgzgd78T06XKy0VvlI0sL9PjmNRW39LSgvHjx2PNmjUAgICAABw7dgwpKSl46KGHzBIg66WuXzjr\nm29M+vGmJuDNN4GXXgIGDhRDX+6/vwvjY2ZhUsIZNGgQxo4dq3ht7Nix+PjjjwEAOp0ORARJkhS1\nD0mS5NsjnU6HxsZG1NXVKWo5kiRBd2WAlE6nM+q1MhgMqK6uVpQpKipSlJEkST53I61LjF6Llxu1\nsN/+FvjqK/H8ttvEClcd9J//iJHC5eVix4TVq9uf+sS6GVMafObOnWvUaLxkyRK699575eP2Go23\nX1msqCONxuXl5WRnZ6doNN6zZ4+i0Xj37t1GjcYbN24kd3d3amxsbDN+bjTuZtLSiF56SSyg1QHV\n1URPPikmT48fT3T4sJnjY13OpIRTVFREjo6O9Mc//pFOnTpFW7dupX79+tG2bdvkMsnJyaRWqyk9\nPZ1KSkooOjqaRo0apegWj4+Pp+HDh1NWVhYVFxdTcHCwUbd4ZGQkBQUFUWFhIR08eJB8fHxo3rx5\n8nmDwUDjxo2jiIgIOnr0KGVkZJBGo6EVK1a0Gz8nnJ6ppYXoww+JNBoiV1eidevaXhWCdX8mr4fz\n+eefk7+/P/Xt25fuuOMO2rRpk1GZVatW0aBBg6hv374UFhZGJ0+eVJy/fPkyJSYm0oABA6hfv34U\nExNDkiQpytTU1FBcXBy5urqSu7s7LViwgOrr6xVlKisrKSoqilxcXEij0VBSUhIZDIZ2Y+eE0/Oc\nOkUUFiZqNQ88QPTf/1o7InYreC4V65YaG8VI4TVrAK0W2LBB9Eaxno1nlLBu58ABYOFCsYbVs88C\nK1f2jtU1GU/etHYo7BrV1WI1hylTxPo0hw4BycmcbHoTvqViVkcEbN0qNpdrbATWrgWeeAKws8k/\nh70b/y9lXc9gAN55R4zKa13btx0nTwLTpwMPPQSEhIg96BYu5GTTW/H/Vlvx5z+LlefGjhUTjczp\nySfFY/VqYMIEsSjNdRoaRIOwvz9w+jSwezewbRsvjtfb8S2VLSgoACZOvHo8dKjYrsBc1Grl1rTv\nvgvMny8f7t8v8tG33wLPPw+sWGG8BRPrnbiGYwt++EF5rNeL7UfMxcenzeOqKuDxx8UuLZ6ewOHD\nYm86Tja2wyYTjs31Uk2dCowcefX48cfN20iSmip2ybzrLmDDBtDkX2LLFrF8xI4donknJwfw8zNf\nCKx74lsqW/Hjj2K3Nw8P4De/6YINmDrm669FI3B2NhAXJ3a4vW5VEWZDeOCfrRgwQNRsLOTyZdG9\n/dprwO23A3v3it4oZts44bAu9+WXolZz5gywdCmwbBnQt6+1o2LdgU224bBbsHOnWPHKxQV4/XXF\nqfPngUceEeNpdDqxW8KaNZxs2FXchtMbNDeLPmaNRrTRmEtTk3j/a7frKS0F3emHf/5TdHEDwJ/+\nBDz6KA/eY8Zs8leiV/VS1dcDv/yl6AIaOhTIyDDftS5fViYbAOWHLmHqVDHM5te/FiOFzd0Jxnou\nruH0dG+/DSQkXD2+806xMLm5/P73wLvv4hKc8OqQDXj93KMYPlyFlBRg2jTzXZb1Dtxo3NNd//fi\nVv5+SJJIYPb2QGJi27dn//gHvhjxJOL/NgZnz7tg2TIV/vAHwMmp85dltoNrOD3dxYuilbawUAzZ\n/de/OrdS1aVLQECAmE0JiOfFxYpNuCVJzOj+8EMxljAlBfD17ZqPwWwD13B6un79gNxckSg0GjHe\npjNOnLiabADRxfTdd8Dw4WhpATZtApKSROVn82bg4YctNnaQ9SLctNcbODiIWeCdTTaAGJ137V4r\nAwYAGg2OHRMLYj3xBDBrlshLjzzCyYZ1jk0mnF7VS9VVPD2Bzz4T2WXaNPy8Yw+WveKMwEAx6TIr\nC3jvPVGMsU4zZcX11atXk0qlUjzGjh2rKPPiiy/KOzaEhoa2uWNDQkKCvGPD7NmzjXZsqK6uprlz\n58o7NsyfP58uXryoKFNZWUkzZswgZ2dn0mq19Pzzz99wxwYi3rWho3bvJvL2JurTh+jll4kuX7Z2\nRKy3MLmG4+fnB0mSoNfrodfrcfDgQflccnIy1q1bh3feeQeFhYVwcXFBeHg4Ghsb5TJLlizB559/\njn//+9/IycnB999/j9mzZyuuMXfuXJSXlyMzMxOff/45cnJy8OSTT8rnW1paMGPGDDQ3NyM/Px9b\ntmzB5s2bsXLlyk6kXNbqhx+AOXPERG9vb6CkBHjxRaBPH2tHxnoNU7LT6tWrKTAwsN3z7e26mZaW\nJh/fbNfNsrIyUqlUil03MzIyFLtu7tq1y2jXzZSUFHJ3d6empqZ24+MaTtsMBqK33yZycyPy9CT6\n4AOx+RxjXc3kGs7JkycxZMgQjBw5EvPmzcPZKyvHVVRUQK/XIyQkRC7r6uqKCRMmIC8vDwBQXFyM\n5uZmRRlfX18MGzZMLpOfnw8PDw95L3IACA0NhUqlQsGVpTHz8/Ph7+8Pz2saFMLDw1FbW4vjx4+b\n+pFsWkkJcO+9QHw88MADYjmJefO4UZiZh0kJZ+LEidi8eTP27NmDlJQUVFRUYMqUKaivr4der4dK\npYL2usVOtFot9Ho9AECSJDg6OhqNgbm2jF6vh0ajUZy3t7eHWq1WlGnrOq3n2M3V14uZ3HffDdTV\niQWx/vEPsTooY+Zi0jic8PBw+bmfnx/Gjx8PLy8vbN++HWPGjOny4Jh57NolZkNIEvDyy8BzzwGO\njtaOitmCWxr45+bmBh8fH5w6dQpTp04FEUGSJEXtQ5Ik+fZIp9OhsbERdXV1ilqOJEnQXVmuX6fT\n4dy5c4rrGAwGVFdXK8oUFRUpykiSJJ+7mTlz5sDBQfnRY2NjERsb29GP3iN9/z2wZIkYjBwWBmRm\nKlceZczsbqUB6MKFC6RWq2ndunVE1H6j8fbt2+XjmzUal5eXk52dnaLReM+ePYpG4927dxs1Gm/c\nuJHc3d2psbGx3XhttdG4uZlo3ToiV1cijYboww+5UZhZh0kJ57nnnqP9+/fTmTNnKDc3l0JDQ0mj\n0VBVVRURESUnJ5Narab09HQqKSmh6OhoGjVqFDU0NMjvER8fT8OHD6esrCwqLi6m4OBgmjx5suI6\nkZGRFBQURIWFhXTw4EHy8fGhefPmyecNBgONGzeOIiIi6OjRo5SRkUEajYZWrFhxw/htMeEcPkx0\nzz1EANGTTxJVV1s7ImbLTEo4c+bMoSFDhpCTkxPdfvvtFBsbS6dPn1aUWbVqlTzwLywsrM2Bf4mJ\nifLAv5iYGKOBfzU1NRQXFycP/FuwYAHV19crylRWVlJUVBS5uLiQRqOhpKQkHvh3jQsXiJ59lsje\nnsjPjyg319oRMUbEs8V7gkuXTFqnc+dOYNEiMSVh1Soxw/u228wYH2MdZJNzqXqM2loxt8nZWazo\nV1Fxw+LffQf89rfAzJlX1+FaupSTDes+bDLh9JjJm2+8ARw4IJ5//bVYH6INBgPwt7+JCeN5eUBa\nmuj6HjHCgrEy1gE2uR5Oampqz7ilqqu78TGAr74S+3QfOiRGC7/6KuDubqH4GDORTdZweoyFC6+u\ncePkJEboXXHhAvDMM8D48WIzhbw8YP16Tjase7PJGk6PMWYMcPy4qL6MGSOmcAP45BNg8WKguhpI\nTgaefprbaVjPwL1UPUhlpUg06elAVJSo0Xh5WTsqxjqOb6msSZJEo3B19Q2LNTcDb74J3HGHWNf8\no49E1zcnG9bT2GTC6Ra9VPn5gI+P6PYeOxb45ps2ixUVAffcI5pvHn8cKC8HZs/m5SNYz8S3VNYy\na5ZojGkVHw9s2CAf1tUBy5eL26aAAGDjRtFAzFhPxo3G1nL9znFXRhITAR9/DDz1lBj398Ybot3G\ngf9PsV7AJm+puoVXXpF7neDvDyxdijNngPvvB2JigF/8AigrE13fnGxYb8G/ytYycqTYeK66Gk1u\nnvjLX1VYvVrsrrtjB/Cb31g7QMa6HtdwrMneHvnfDsQv7hH7cz/xhGgU5mTDeiubTDjdoZfqp5/E\nMp/BwWLQXmEh8NZbQP/+VguJMbPjXioLIxJLfD79NHDxopj7tGiR2LObsd7OJms41lJRAcyYAfzu\nd8CkSeL26amnONkw28EJxwKamsScpzvvFFOjPv1UdH0PHWrtyBizLO6lMrP//EcsH1FeLm6jXnoJ\n6NfP2lExZh1cwzGTmhqRaO69VyzYV1wsBvFxsmG2zCYTjjl7qYiAbdvEahLbtgHr1olazl13dfml\nGOt5bmUF9tdee41UKhU988wzitdffPFFeeeG0NDQNnduSEhIkHdumD17ttHODdXV1TR37lx554b5\n8+fTxYsXFWUqKytpxowZ5OzsTFqtlp5//vkb7txg7l0bTp0imj5dbMnywANE//2vWS7DWI/V6YRT\nWFhI3t7edNdddykSztq1a8nDw4N27txJpaWlFB0dTSNGjFDsTbVw4ULy8vKi7OxsOnToEE2aNMlo\nb6qIiAgKDAykoqIiys3NpdGjR1NcXJx83mAwkJ+fH4WFhVFJSQllZGTQwIEDafny5e3GbK6E09BA\n9OqrRE5ORF5eRJ991qVvz1iv0amEc+HCBfLx8aHMzEyaOnWqIuG0t/tmWlqafHyz3TfLyspIpVIp\ndt/MyMhQ7L65a9cuo903U1JSyN3dnZqamtqM2xwJJyeH6I47xP5PSUlE11XCGGPX6FQbzqJFi3D/\n/fdj2rRpitcrKiqg1+sREhIiv+bq6ooJEyYgLy8PAFBcXIzm5mZFGV9fXwwbNkwuk5+fDw8PD3lP\ncgAIDQ2FSqVCQUGBXMbf3x+enp5ymfDwcNTW1uL48eOd+Vgmqa4Gfv97sZxN//5iFdDkZMDFBcDl\ny8CLLwJz54r+b8YYgE50i6empuLIkSMoLi42OqfX66FSqaDVahWva7Va6PV6AIAkSXB0dDQa6Xtt\nGb1eD41Gozhvb28PtVqtKNPWdVrPBQQEmPrROoQI2LpVbC7X2Ai8/baYA2V3bepOTAQ2bRLPU1OB\nzEzgvvvMEg9jPYlJNZzvvvsOS5YswdatW3GbDa7a/c03wPTpwEMPAdOmibE1Cxdel2wAIDf36nMi\n0U3FGDOthvPVV1/h/PnzuPvuu0FXpmAZDAbk5ORg3bp1OHHiBIgIkiQpah+SJMm3RzqdDo2Njair\nq1PUciRJgk6nk8ucO3dOcW2DwYDq6mpFmaKiIkUZSZLkczcyZ84cOFy3yExsbCxiY2PbLN/QALz+\nupj3NHgwsHs3EBFxgwtMmgScOCGeq1TimDFmWrf4xYsX6fjx44rHPffcQw8//DCVlZURUfuNxtu3\nb5ePb9ZoXF5eTnZ2dopG4z179igajXfv3m3UaLxx40Zyd3enxsbGNuPvTKNxdjaRry+RgwPRCy8Q\n1dd34Id+/lkUfvBBon/9q8PXYqy3u6VxOERk1EuVnJxMarWa0tPTqaSkhKKjo2nUqFGKbvH4+Hga\nPnw4ZWVlUXFxMQUHBxt1i0dGRlJQUBAVFhbSwYMHycfHh+bNmyefNxgMNG7cOIqIiKCjR49SRkYG\naTQaWrFiRbuxmpJwzp8nevRRMaYmOJiotNSUb4Ux1pZbTjj33Xef0cC/VatWyQP/wsLC2hz4l5iY\nKA/8i4mJMRr4V1NTQ3FxcfLAvwULFlD9ddWLyspKioqKIhcXF9JoNJSUlHTLA/9aWog2byYaMIDI\n3Z3onXeIbvCWjDET8Ho41/j6a9EInJ0NxMWJuU/XdYQxxm6BTc6lut7ly8Dq1cC4ccDZs8DevcD/\n/i8nG8a6mk0uT9HaSxUbGwutNhYLFwJnzgBLlwLLlsk7tjDGuphN1nBSU1OxaVM6du+ORUgIoNMB\nR44Aa9Z0INnk5QEjRojhxStWWCRexnoLm2zD+fvfa7FqlSuIgD//GXj00TYG77Vn+HDg//7v6vGX\nX/IoYsY6yKZqOK1j8RYvBqKixPHjj5uQbADgugGJRseMsXbZTMIhEivwAUB6OvD++8B107U6JjHx\n6vPRo4Hw8C6JjzFbYFO3VCUldQgI6IJtYr74AqiqAiIjAXf3rguQsV7OphJOaxtOZGSk3EvV3vwp\nxljXs8mEY82N8BizZTbThsMYsz5OOIwxi+GEwxizGE44jDGLscmEY86N8Bhj7eNeKsaYxdhkDYcx\nZh2ccK6l1wPPPAMkJAAnT1o7GsZ6HZtcD6dNBgMQEgKUlYnjHTvEPjA8dYGxLsM1nFaSdDXZAKK2\n0zq9nDHWJWwy4bTZSzVwIODldfXY3R0YNcrywTHWi5mUcFJSUhAQEAA3Nze4ubkhODgYGRkZijIr\nV67E4MFqhb/vAAAJs0lEQVSD4ezsjOnTp+PUqVOK8w0NDVi0aBE8PT3Rv39/xMTEGG16V1NTg7i4\nOLi5ucHDwwMLFixAfX29oszZs2cRFRUFFxcX6HQ6JCUloaWlpUOfIzU1Fenp6cqJm7fdJmaBz5kD\nzJolFja+Zt9yxlgXMGWLh88++4x2795Np06dopMnT9Ly5cvJ0dFR3gRv7dq15OHhQTt37qTS0lKK\njo6mESNGKPakWrhwIXl5eVF2djYdOnSIJk2aZLQnVUREBAUGBlJRURHl5ubS6NGjKS4uTj5vMBjI\nz8+PwsLCqKSkhDIyMmjgwIG0fPnyG8bfmY3wGGNd55b3pVKr1fTee+8RUfu7bqalpcnHN9t1s6ys\njFQqlWLXzYyMDMWum7t27TLadTMlJYXc3d2pqamp3Vg54TBmXZ1uw2lpaUFqaip+/vlnBAcHo6Ki\nAnq9HiEhIXIZV1dXTJgwAXl5eQCA4uJiNDc3K8r4+vpi2LBhcpn8/Hx4eHjIe5EDQGhoKFQqFQoK\nCuQy/v7+8Lzmlic8PBy1tbU4fvx4Zz8SY8zMTE44x44dQ//+/dGnTx8kJCRgx44d8PX1hV6vh0ql\ngva6zZy0Wi30ej0AQJIkODo6Go3yvbaMXq+H5rq1P+3t7aFWqxVl2rpO6znGWPdk8jicMWPG4OjR\no6itrcVHH32Ehx9+GDk5OeaIzWxa96W6Fq/+x5j5mZxwHBwcMGLECABAYGAgCgsL8de//hVJSUkg\nIkiSpKh9SJIk3x7pdDo0Njairq5OUcuRJAk6nU4uc32vlcFgQHV1taJMUVGRoowkSfK5Np04ASQn\nAwBS//IXuPr4mPrRGWO36JbH4bS0tKChoQHe3t7Q6XTIzMyUz9XV1aGgoADBwcEAgKCgIDg4OCjK\nfP3116isrMSkSZMAAJMmTcJPP/2Ew4cPy2UyMzNBRJgwYYJcprS0FFVVVXKZvXv3ws3NDXfccYdx\nkNXVwK9+BWzeLI5nzgQ62IXOGOtCprQwv/DCC5STk0Nnzpyh0tJS+sMf/kD29vaUmZlJRETJycmk\nVqspPT2dSkpKKDo6mkaNGqXoFo+Pj6fhw4dTVlYWFRcXU3BwsFG3eGRkJAUFBVFhYSEdPHiQfHx8\naN68efJ5g8FA48aNo4iICDp69ChlZGSQRqOhFStWtB34wYNEANUCopcKINLrTfnojLEuYFLCmT9/\nPnl7e5OTkxNptVqaPn26nGxarVq1igYNGkR9+/alsLAwOnnypOL85cuXKTExkQYMGED9+vWjmJgY\nkiRJUaampobi4uLI1dWV3N3dacGCBVRfX68oU1lZSVFRUeTi4kIajYaSkpLIYDC0Hfi5c0Tu7lcT\njpcXUXOzKR+dMdYFbGc9nOJi1L3yCtw+/RS1R47ANSDA2hExZnNsJ+GAF+BizNpscvImY8w6bKqG\nQ0S4cOEC+vfvD5VKZe1wGLM5NpVwGGPWxbdUjDGL4YTDGLMYTjiMMYvhhMMYsxibSzg9dbfNnhh3\nT4wZ4LjNiRNOD9ET4+6JMQMctznZXMJhjFkPJ5w2dOQvhSXLdFRPjLsrr9UT4+7NvyNt4YTTBv5l\n6nn/cLv6vSx1rd78O9KWXr3Vb+tUhms1Nzejrq7uhj/X3cp0x5i6W5nuGJOtff6OTBnq1VMbWmeH\nM8bMryOrMPTqhNNWDYcxZh42X8NhjHUv3GjMGLMYTjiMMYvhhMMYsxhOOIwxi+GEwxizGJtJOOvX\nr4e3tzf69u2LiRMnGm0V3N289NJLsLOzUzza3FXUyg4cOICZM2diyJAhsLOzQ3p6ulGZlStXYvDg\nwXB2dsb06dNx6tQpK0SqdLO4H3vsMaPvf8aMGVaKVnjttdcwfvx4uLq6QqvVYtasWfjmm2+MynXH\n77uVTSSctLQ0PPvss3jppZdw+PBhBAQEIDw8XLFVcHfk5+cHSZKg1+uh1+tx8OBBa4dkpL6+Hnfd\ndRc2bNjQ5hiM5ORkrFu3Du+88w4KCwvh4uKC8PBwNDY2WiHaq24WNwBERkYqvn9rz8Y+cOAAFi9e\njIKCAuzbtw9NTU0ICwvDpUuX5DLd9fuWWXzrPSuYMGECPfXUU/JxS0sLDRkyhJKTk60Y1Y2tXr2a\nAgMDrR2GSVQqFX366aeK1wYNGkRvvvmmfFxbW0tOTk6UlpZm6fDa1Vbcjz76KM2aNctKEXXM+fPn\nSaVS0YEDB+TXuvv33etrOE1NTfjqq68QEhIiv6ZSqRAaGoq8vDwrRnZzJ0+exJAhQzBy5EjMmzcP\nZ8+etXZIJqmoqIBer1d8966urpgwYUK3/+4BIDs7G1qtFmPGjEFCQgKqq6utHZLCTz/9BJVKBbVa\nDaBnfN+9PuFUVVXBYDBAq9UqXtdqtdDr9VaK6uYmTpyIzZs3Y8+ePUhJSUFFRQWmTJmC+vp6a4fW\nYXq9HiqVqsd994C4nXr//ffx5Zdf4vXXX8f+/fsxY8YMUDcZmE9EWLJkCSZPniy37fWE77tXzxbv\nycLDw+Xnfn5+GD9+PLy8vLB9+3Y89thjVozMNjz44IPy8zvvvBP+/v4YOXIksrOzcd9991kxMiEh\nIQFlZWXIzc21digm6fU1HE9PT9jb20OSJMXrkiRBp9NZKSrTubm5wcfHp1v1ONyMTqcDEfX47x4A\nvL294enp2S2+/8TEROzatQvZ2dkYNGiQ/HpP+L57fcK57bbbEBQUhMzMTPk1IkJmZiaCg4OtGJlp\nLl68iG+//VbxC9bdeXt7Q6fTKb77uro6FBQU9KjvHgC+++47/Pjjj1b//hMTE/Hpp58iKysLw4YN\nU5zrEd+3VZusLSQtLY369u1LW7ZsofLycnriiSdIrVbTuXPnrB1au5577jnav38/nTlzhnJzcyk0\nNJQ0Gg1VVVVZOzSFixcv0pEjR+jw4cOkUqnorbfeoiNHjlBlZSURESUnJ5Narab09HQqKSmh6Oho\nGjVqFDU0NHTbuC9evEjPP/885efn05kzZ2jfvn0UFBREY8aMocbGRqvFHB8fT+7u7pSTk0N6vV5+\nXLp0SS7TXb/vVjaRcIiI1q9fT15eXuTk5EQTJ06koqIia4d0Q3PmzKEhQ4aQk5MT3X777RQbG0un\nT5+2dlhGsrOzSaVSkZ2dneLx2GOPyWVWrVpFgwYNor59+1JYWBidPHnSihELN4r70qVLFB4eTlqt\nlvr06UPe3t60cOFCq/+BaiteOzs72rJli6Jcd/y+W/F6OIwxi+n1bTiMse6DEw5jzGI44TDGLIYT\nDmPMYjjhMMYshhMOY8xiOOEwxiyGEw5jzGI44TDGLIYTDmPMYjjhMMYs5v8BZ7he/wn+QvoAAAAA\nSUVORK5CYII=\n", "text/plain": [ "Graphics object consisting of 2 graphics primitives" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "gtPlot = list_plot(GtList, rgbcolor='red', figsize=3)\n", "gtFunc(x) = GtModel.subs(GtFit)\n", "calGtPlot = plot(gtFunc, [x, 0, 19], figsize=3)\n", "(calGtPlot + gtPlot).show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

1967年から1989年までの予測

\n", "\t

\n", "\t\tちょっと乱暴な遊びですが、\n", "\t\t政府支出関数gtFuncを使って、1967年から1989年まで政府支出を求め、前年度の消費Ct1List\n", "\t\tの値とgtFuncの値から消費、投資、所得を順番に計算してみましょう。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# 予測\n", "predCt1List =[ Ct1List[0]] + range(23)\n", "predCtList = [calCtList[0]]\n", "predItList = [calItList[0]]\n", "predYtList = [calYtList[0]]\n", "for i in range(1, 24):\n", " predCt1List[i] = ctFunc(predCt1List[i-1], gtFunc(i-1))\n", " predCtList.append(ctFunc(predCt1List[i], gtFunc(i)))\n", " predItList.append(itFunc(predCt1List[i], gtFunc(i)))\n", " predYtList.append(ytFunc(predCt1List[i], gtFunc(i)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\t

予測結果

\n", "\t

\n", "\t\t政府支出の乱暴な予測にも関わらず、予測値(緑)は先のシミュレーション結果と\n", "\t\tほぼ同じ傾向を示しています。\n", "\t

\n", "\t

\n", "\t\tこのようにシミュレーションは、未来の予測にも活用できることが理解して頂けたでしょうか。\n", "\t

\n", "" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false, "scrolled": true }, "outputs": [ { "data": { "text/html": [ "
政府支出消費
投資所得
" ], "text/plain": [ "<__main__.Table2Html instance at 0x7fa32655d050>" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# 計算値(青)、予測値(緑)と実測値(赤)\n", "predCtPlot = list_plot(zip(range(24), predCtList), plotjoined=true, rgbcolor ='green', linestyle='dashed', figsize=3)\n", "predItPlot = list_plot(zip(range(24), predItList), plotjoined=true, rgbcolor ='green', linestyle='dashed', figsize=3)\n", "predYtPlot = list_plot(zip(range(24), predYtList), plotjoined=true, rgbcolor ='green', linestyle='dashed', figsize=3)\n", "calGtPlot = plot(gtFunc, [x, 0, 23], figsize=3)\n", "# プロット\n", "ctPlt = (calCtPlot + predCtPlot + ctPlot + valCtPlot)\n", "itPlt = (calItPlot + predItPlot + itPlot + valItPlot)\n", "ytPlt = (calYtPlot + predYtPlot + ytPlot + valYtPlot)\n", "Table2Html([[\"政府支出\", \"消費\"], [_to_png(calGtPlot + gtPlot), _to_png(ctPlt)], \n", " [\"投資\", \"所得\"], [_to_png(itPlt), _to_png(ytPlt)]])" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "SageMath 7.2", "language": "", "name": "sagemath" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 2 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", "version": "2.7.10" } }, "nbformat": 4, "nbformat_minor": 0 }