Mathematica tutorial for young economic theorists:
応用経済理論屋のための Mathematica ガイド
version 1.4
大城 淳
∗2017 年 12 月 18 日
概要
本稿は経済学，特に応用ミクロ経済学(マクロでも可)で論文を書こうとする人向
けに，Mathematica⃝Rの使い方をメモしたものである．およそ，Mathematica歴1 週間以上・プログラミング経験ゼロ近傍の文系読者を想定している．
1 はじめに
経済学研究者にとってほどよい Mathematica のテキストはなかなか見当たらない． Mathematicaはネット上のリソースも少ない．RやPythonと比べると，Mathematica は遠からず滅ぶのではないかとも危惧する．本稿は，経済学界隈でMathematicaのさら なる普及を願って，私のつたないテクニックを書き留めたものである．誰しもソフトウェ アの習熟にはなるべく時間を割きたくはないもので，その助けになれば幸いである．よく 使う操作・関数を眺めながら，モデルを解く，グラフを書く，統計データをいじる，など 分野を問わず経済学全般で参考にできそうな話を，包括的にはカバーできていないが，ま とめてある．私はプログラミングのプロでもなければ研究歴も浅いのだが，いくつか論文 を書く中で学んだことを書き留めておきたい．
Mathematica 初心者の院生は，指導教員や先輩などのnotebook を見せてもらおう．
自分で編み出すよりも効率的に技能を伸ばすことができるはずである．シンタックスは異
なるが，Matlab のコードも参考になる．Matlab の場合何よりネット上に出回っている
∗Department of Law and Economics, Okinawa University. j-oshiro@okinawa-u.ac.jp.
ことが多い．もちろん本稿に出てくるコードを実行しつつ解読していくのも勉強になるは ずである．
初歩的な疑問は，ヘルプかGoogleを頼ればたいてい解決する．Mathematicaのヘルプ はよくできており，ヘルプの例題を根気強く追いかけていれば自然と上達する．ヘルプと
Googleをまめに参照することは当然であろうから，本稿では関数の使い方について細か
い説明は省く．ヘルプよりも詳しいものは，Wolfram Mathematica Tutorial Collection https://www.wolfram.com/learningcenter/tutorialcollection/
を参照するとよい(が実際に参照する機会は少ないかも)． 参考書を買うなら，どんぴしゃりのものはないが強いて言えば
(1) Sal Mangano (2011) 『Mathematicaクックブック』 オライリージャパン (2) 谷口義治・永友清和 (2012) 『線形代数とMathematica』 牧野書店
(3) 上坂吉則 (1997) 『Mathematica数値数式プログラミング』 牧野書店
の3冊が有用である．理系向けのものにはもっとごついものもあるが，数学で勝負したい 人以外には過剰だと思う．
本稿と趣旨が近いもので，Nasser M. Abbasi氏のチートシートは便利である． http://12000.org/my notes/faq/mma notes/MMA.htm
本稿はversion 11.0 (Windows版)に基づいて執筆している．多少古いバージョンでも 問題ないかと思う．
1.1 他のソフトウェアと比べて
Mathematicaは応用理論で戦うべく代数の計算をせっせとこなす必要があるような人
にはうってつけのソフトである．高額だが理論屋なら優先して購入すべきだと思う． Mathematicaと競合するソフトに，Maple, MathStudioやフリーウェアのMaximaが ある．残念ながらこれらは身近にユーザーがほぼおらず，学ぶメリットがない．Maxima は素人が触るとすぐクラッシュするのでおすすめしにくい．
実証分析にはR, Stata, Pythonなどを使おう．Mathematicaでもそれなりにデータは 扱えるが，ちょっとした検定やちょっと手の込んだ標準偏差を計算するにも一苦労する．
Mathematica自体の計算速度が遅いこともさることながら，何より計算に必要なコード
を書く時間が無駄である．Mathematicaで生成したデータは早くcsv にExport して適
当な統計用ソフトに移行した方がよい．それでもなお，私はMathematicaによる統計的 分析をほんのり経験したので，第5章では経済学で使いそうな統計分析やそれに伴う行列 の操作について少しだけ触れる．
2 入門編
2.1 覚えておくべき操作
2.1.1 背景色
真っ白い背景色は目が疲れる．「書式 (R)→スタイルシート(Y)→ReverseColor(色の 反転)」で，背景を黒くしておくと目に優しい．真っ黒だと落ち着かない場合，「編集(E)→ 環境設定(S)→詳細→オプションインスペクタを開く→書式設定」から，notebookの
背景色(Background) を目に優しい色(灰色など)にしてもよい．ブルーライトカットを
する他のソフトウェアを併用して目をいたわるべし．
2.1.2 中断
Windowsの場合，「Altキー + ピリオド.」のショートカットを覚えよう．いつ終わる
かわからない計算を中断させたい場合にはすかさず使おう．それでも無理そうなときは，
「評価(V)→カーネルを終了(Q)→Local(ローカル)」でいったん振り出しに戻す．
Mathematicaではカーネルと呼ばれるものがあらゆる計算を担っている．ユーザーの
インプットを解釈してアウトプットするエンジンだ．アウトプットはnotebookの左側に Out[(数字)]とナンバリングされる．
2.1.3 %
直前の結果を参照したいときは %を使う．パーセントを重ねていけばもっと前までた どることができ，たとえば%%%とすれば3つ前の結果が得られる．%nとすれば，Out[n] を参照することができる．変数名を定義するほどでもない計算をするときによく使う．
2.1.4 前置修飾・後置修飾
多くの場合，関数は重ねて使う．たとえば
Flatten[Simplify[Solve[x == y, z]]]
という具合だ．上のコードは，次のコードと等価であると知っておくと，すっきりしたプ ログラムを書くことができる:
Flatten@Simplify@Solve[x == y, z]
または
Solve[x == y, z] // Simplify // Flatten
要は，@が前置修飾に使え，//が後置修飾に使えるのである．入れ子構造を解消して読 みやすくなるばかりか，コードに流れができ理解しやすくなる．
2.1.5 コメント
プログラムにコメントを挟みたいときは (* ほげほげ *)とする．Windows の場合
「Altキー + スラッシュ/」というショートカットを覚えておけばコメントアウト・解除 がとても楽になる．
ある程度コメントを書き込んでおかないと，共著者はおろか，半年後の自分さえ解読に 苦労する羽目になる．コメントは長期的に活きるのでこまめに使おう．
クオーテーションを使うと，評価の合間に文章を挟むことができ，読みやすいnotebook になる．たとえば
"the number of regions" iI = 47
"the length of time" tT = 10
"parameters following Hayashi-Prescott"
para = {theta -> 0.362, delta -> 0.089, beta -> 0.976, alpha -> 1.373, tau -> 0.480}
という具合だ．
長い文章を書き込んでおきたい場合は，notebook右端のグループブラケットをクリッ クしてから，スタイルをText (Altキー + 7)にして，文章だけのセルを確保する．
日本語は，フォントをいじくってるマシンではうまく表示されなかったり文字化けし たりもする．苦痛でない範囲で，コメントは英語を使用するくせをつけておいたほうが よい．
2.1.6 出力の省略
計算結果をわざわざ出力するまでもない場合は，セミコロン;を使って抑制する． 出力はするけど全部はいらない，というときはShortを使うと短く表示できる．
2.1.7 変数名の定義
Mathematicaでは大文字・小文字が区別される．標準装備の関数は，先頭が大文字に
なっている．それとの競合を避けるため，自分でシンボルを定義する場合は小文字でス タートするとよい．たとえば上の例では，地域数をiI = 47としているが，これはI だ けだと虚数単位として認識されてしまうためだ．
変数名は好きなあの子の名前でもなんでもよいのだが，論文中の表記と合わせて，それ
なりにmake senseなものにしておいたほうが無難だ．少々長くても，Tabキーで入力補
完してくれるのだからあまり気にならないはず．
変数名には他にもいくつか注意がある: 数字で初めてはいけない．x2は変数名として
OKだが，2xは2 × xと解釈されてしまう．ギリシャ文字は円周率πそのものを除き利
用することができる(π2などはOK)．ESCキーで挟んで入力すると早い． 添え字は定義に利用することができる．たとえば下付の添え字を使い
xn := N[(1 + 1/n)^n]
とすると，x1000= 2.71692などと出力されるようになる．
2.1.8 [[…]], Flatten
Mathematica では主にリストと呼ばれるものを扱う．リストは波括弧でくくられ，
{expr}，という形をしている．
リストの要素にアクセスするには，二重角括弧を使う．たとえば，x を長さ3のベク トル{a,b,c}とする．このとき x[[3]] とすると，c を得る．またたとえば行列A のi 行j 列成分にアクセスしたい場合，A[[i,j]], と書けばよい．i 行だけを取り出したい 場合，A[[i, All]] と書く．i 行のうち，j = 3, 4, . . . , 10, の部分を抽出したければ， A[[i, 3;;10]], と書く．Pythonと違い，番号は1からスタートする．リストの階層が 深くなっても，同様にA[[i,j,k,h]], などと書いていく．詳しくは関数Partのヘルプ を参照せよ．
個人的に，Mathematicaでよくあるエラーの原因は，スカラーとベクトルが混在してし まうことである．たとえばxを使っているつもりが{x}を使っていた，などだ．後者を前者 に直すには{x}[[1]]やFirst@{x}とすればよい(ただし，xが空集合だとエラーが出る)． 具体的にたとえば条件分岐を使ったプログラムを書く際に，If[{1} > 0, t, f]のよう に比較できず真偽を判断できないような書き方だと評価されずうまくいかない．なお，不 等号はベクトルの比較をそのままではしてくれないので，{1}>{0}や{1,2}>{0,0}などと
書いてもTrueとはならず評価されない．一般的なn-ary relation(二項関係≿など)を定 義したければ自分で関数を書く．
波括弧を外す際にはFlattenが活躍する．たとえばFlatten[{{{{{x, {y, z}}}}}}] とすれば{x,y,z}を得る．
2.1.9 小数点
「0」と「0.」と「0.0」はMathematicaでは意味が異なる(他のソフトウェアでも同様)．
「0」のように単に整数を書いた場合は厳密な値として扱われる．ところが「0.」のようにピ リオドがついた場合，実数の近似値として扱われ，整数ではなくなる．また「0.0」は1桁 の精度を持つ近似値という意味になる．たとえば，整数かどうか判定する関数IntegerQ を使って
IntegerQ[1] IntegerQ[1.]
と書くと，前者はTrueが返ってくるが後者は Falseが返ってくる．IfやWhileなどで 条件分岐を考えるときにこうした違いでひっかかることがある．たとえば以下の2つを比 べよう．
If[1∈ Integers, t, f] If[1.0 ∈ Integers, t, f]
前者は1は整数なので tを返す．ところが後者は評価されず，入力したそのままのもの が返ってくる(f とは返ってこない! If[IntegerQ[1.0], t, f]と書けば f と出力さ れる)．
3 定番の関数
Mathematica を経済学で使う主な目的は，(1) 式を整理する，(2) 方程式を解く，(3) 図を書く，の3パターンに分けられる．それぞれ代表的な関数を簡単に紹介する．
3.1 式を整理する
3.1.1 Simplify
Simplifyはその名の通り，式を整理してくれる．使いこなしさえすればとても便利だ．
FullSimplifyはもう少し時間をかけていろいろ試してくれる．ややこしい関数が出てき たらFullSimplifyを試してみるべきだ．
しかしこれらは自分が思った通りに変形してくれるとは限らない．特に，指数が混じっ ているときには要注意だ．変数の定義域を指定するとうまくいくことがある．定義域を指 定するには，
Simplify[z, x > 0 && y > 0]
のように，コンマの後に条件を加えていけばよい．ここではxとyが正，という条件を書 いている．もっと複雑な条件を課すことも可能だ．実践上は，ある程度制約を課さないと 思い通りにSimplifyが機能しないことが多い．経済学で扱う変数の多くは正の実数なわ けで，虚数の可能性までしっかり加味してくれるのはありがた迷惑でもある．
たとえば簡単な効用最大化問題を解く際に出てくる一階条件を見てみよう． D[x^a y^(1 - a) + lambda (M - x - p y), x]
Simplify@Solve[% == 0, x] この結果は，
{{x →( y
(-1 + a)lambda
a
)-1 + a1
}}
という具合になる．(ya−1)1/(a−1) = yのように指数の部分を簡約化するにはもう少し仮 定が必要なのだ．こういうときには，
D[x^a y^(1 - a) + lambda (M - x - p y), x] Simplify[Solve[% == 0, x] , a < 1 && y > 0] と仮定をおいておけば，無事
{{x → y( lambda a
)-1 + a1 }} となる．
仮定を置いて式を整理するには，Refineを使用する手もある．たとえば，上の2行目 で使ったSimplifyをRefineに置き換えても同じような出力を得る．
Simplify に制約をつけなくても，PowerExpandで無理矢理展開すればうまくいくこ とがある．PowerExpandは中身が全部正の実数だと仮定して展開してくれる．もちろん． 実際に正の実数が含まれているかどうかわからないときに PowerExpandを使うのは危 ない．
丁寧かつ発見的に上の効用最大化問題を解くなら，たとえば次のようなコードになるだ ろう．
util[x_, y_] := x^a y^(1 - a) budget[x_, y_] := M - x - p y
lag[lambda_, x_, y_] := util[x, y] + lambda budget[x, y] (* dL/dx = 0 *)
foc1 = Flatten@Simplify[
Quiet@Solve[D[lag[lambda, x, y], x] == 0, x], y > 0 && a < 1 && M > p y]
(* combining this with dL/d lambda = 0, or budget = 0 *) Solve[D[lag[lambda, x, y], lambda] == 0, x][[1, 1, 2]] foc2 = Flatten@Simplify@Quiet@Solve[% == x /. foc1, lambda] (* solve for y *)
Simplify[Solve[D[lag[lambda, x, y], y] == 0, y][[1, 1, 2]] /. foc1 /. foc2, 1 > a > 0 && y > 0 && M > p y]
Flatten@Simplify[Quiet@Solve[% == y, y], 1 > a > 0 && y > 0]
solution = Simplify[{x, y, lambda} /. foc1 /. foc2 /. % , 1 > a > 0 ]
Quietは余計なエラーメッセージを黙らせるものであり，ここでは不要．エラーをはか
なくなるまでループを作って計算させるときなどに使うことがある．また，効用関数や 予算制約やラグランジュ関数をx, y, λの関数であることがわかるように明示しているが， ここではその必要はない．
経済学の問題を解くときには，ある程度手計算によるフォローが必要になるのはやむを 得ない．むしろ，よくあるCobb-DouglasやCESの最適化ぐらいなら手計算したほうが 早いと思う．上のコードから冗長な部分をなくし，もう少し人間が解くときのようにする と次のようになるだろう:
D[lag[lambda, x, y], {{x, y}}] // Simplify
%[[1]] p - %[[2]] // Simplify (* eliminating lambda *) foc1 = Solve[% == 0, x] // Simplify // Flatten
foc2 = Solve[budget[x, y] == 0 /. foc1, y] // Simplify // Flatten foc1 = foc1 /. foc2 // Simplify
solution = Flatten@{foc1, foc2}
上から2行目は，λ = . . . とpλ = . . . を連立して解いているだけ．ここは
% . {1, -1/p} // Simplify (* eliminating lambda *)
などと内積を使って書いてもよい(ドット.は行列のかけ算にも使う)．また，下から2行 目でfoc1を上書きしているが，上書きせず別のシンボルに格納してもよいしそのほうが 混乱は少ないかもしれない．
多少式が汚くても，計算していくうちにきれいになることもある．逐一手計算さながら に式を整理しなくてもよい，と考えると楽だ．SimplifyやFullSimplifyはそれなりに 計算時間を食うので，必要がなければセミコロンで出力を抑制しただけで済ませておけば
よい(第2.1.6章)．特に何千何万回も評価するプログラムの中には入れないほうがよい．
3.1.2 Factor, Expand
Simplifyの補助として，因数分解Factorと，展開Expandを使うことがある．上述
の通りPowerExpandも使用できる場面は少ないが便利である．
思い通りの形に因数分解してくれないことは多いが，式の整理にあまり時間を費やさず 突き進むとよいと思う．モデルの特定化自体を見直したほうが早いこともある．
3.1.3 CoefficientList
複雑に見える多項式を整理していきたいときがある．たとえば
(ちょー長いパラメーターの塊)x2+ (塊’)x + (定数項) = 0,
の各係数部分をもう少し意味のわかるようにしたいことがままある．こういう場合，
CoefficientListを使えば，特定の項の係数だけを取り出してくることができる．連立
方程式を行列で表現したいときに便利である．部分的にSimplifyを使ったり，何らかの パラメータの値をゼロにしてみたりすることで，思いの外わかりやすい形に変形できるこ とがある．
二変数ある場合，CoefficientList[(何らかの多項式), {x,y}]のようにして使うと， 行列が出力される．行で見て，上から1の係数，xの係数, x2の係数，…となる．列で見 て，左から1の係数, yの係数, y2 の係数，になる．たとえば(2,2)要素は，xyの係数で， (3,6)要素はx2y5 の係数，という具合だ．
3.1.4 Limit
極限における値を知りたい場合，Limitを使う．近づく方向はDirectionで決めるこ とができる．
3.1.5 Piecewise, Max, Min, UnitStep
変数の定義域がややこしいときや関数が不連続なときに数値計算を行う場合，普通の ニュートン法では収束する保証はない．場合分けがあるときには，Piecewiseで作成す ることができる．不連続関数の数値計算がうまくいくこともある．
場合分けを処理するには，Max, Min, If, のような関数もうまく組み合わせる．
Piecewiseの代わりに/;を使って書くこともできる．Ifを使って書くよりも可読性 が高くなると思う．たとえば，プロスペクト理論の価値関数を図示してみよう．
v(x) =
{x0.5 if x ≥ 0
−2.5(−x)0.5 otherwise この式は次のように定義することができる:
(*---define value function v---*) v[x_] := x^0.5 /; x >= 0
v[x_] := - 2.5 (-x)^0.5 /; x < 0 (*---illustration---*) Plot[v[x], {x, -100, 100},
PlotRange -> {{-100, 100}, {-20, 20}}, (* limit of both axes (xlim & ylim) *) AxesLabel -> {x, v[x]} (* labels for axes *) ]
-100 -50 50 100x
-20 -10 10 20 v(x)
Piecewiseの特殊ケースとして，x < 0のとき0, x ≥ 0のとき1，となるstep function はUnitStepで書くことができる．符号を判定して1,0,-1,のいずれかを割り当てるSign
や，値域を0,1で切り取るClipも似たようなものである．これらは条件分岐や離散選択 モデルなどで使える．
3.2 方程式を解く
3.2.1 D
微分はDを使う．2階微分は， D[f[x], {x, 2}]
とすればd2f /dx2 を計算してくれる． 交差微分∂2f /∂x∂yは，
D[f[x, y], {x, y}] とする．一方，
D[f[x, y], {{x, y}}]
は{∂f/∂x, ∂f/∂y}のようにベクトルでの微分を求める．コードが似ているので使い分
けに注意．
ヘッシアンを計算するには，
D[f[x, y], {{x, y}, 2}]
とすればよい．行列はぱっと見で判断しにくいので，MatrixFormをかましてコードが正 しいかどうか確認するとよい．二階条件の確認に行列式を評価する必要があればDet を 用いる．
テイラー展開をしたい場合，Seriesを用いる．たとえば Normal@Series[f[k, l], {k, k0, 2}, {l, l0, 2}]
は，関数f (k, l)を(k, l) = (k0, l0)の周りで二次まで近似したものを出力する．Normal を省くと，高次の項o(x − x0)3を含む表現となる．
積分はIntegrateを，数値的に求める際はNIntegrateを使う．フーリエ変換やラプ ラス変換のような積分変換用の関数も用意されている(が使ったことはない)．
3.2.2 Solve, Reduce
読んで字のごとくSolveで方程式を解くことができる．連立方程式も解ける．たとえ ば，∂f /∂x = ∂f /∂y = 0という連立方程式を解くなら，
Solve[D[f, {{x,y}}] == 0, {x,y}]
とする．もちろん，内点解以外を求めたいときにも多用する．
Solveの解は，{{x → 1}}のように，矢印(Rule)付きで出てくる．矢印には興味がな い，という場合，
x /. Solve[f == 0, x]
のように，xに代入してしまうか， Solve[f == 0, x][[1,1,2]]
と矢印の右側部分だけ取り出すかするとよい．ただ，矢印付きのまま扱っても害はないど ころか便利なことも多く，逐一矢印を排除しなくてもよい．
個人的には，SolveはFlattenといっしょに使い，かっこの階層があまり深くならな いようにすることが多い．
Solve と似たものにReduceがある(私はあまり使ったことがない)．解が満たすべき 制約条件を課して解きたいときに使える他，Reduce は特殊解も教えてくれる．出力形 式をSolveのような規則リスト(x → 1. みたいな矢印付きのやつ)に直したいときは， ToRulesを使用する．
3.2.3 NSolve, FindRoot
数値的に解を求めたい場合，Solveが機能しないことがある．そういう場合は，NSolve ないしFindRootを使う．前者は複数(だが少数)の解を知りたいときに，後者は意味の ある解を一つだけ知れば十分なときに使う．具体的な数値計算プロセスについては，4章 で触れる．
複雑なシミュレーションをするときには NSolve では遅すぎるので，FindRoot を使
う．FindRootはニュートン法の類いのアルゴリズムを用いており，結果は初期値に依存
する(NSolve はグレブナー基底を使って計算する)．非凸関数を解析する場合に限らず，
FindRootの初期値は何パターンか試した方がよい．解に十分近いところからスタートし
ているはずなのになかなか解にたどり着いてくれないことはよくある．グラフを書いてみ
て解のありそうな範囲を特定する，複素根の可能性を見当する，など工夫をこらすべき だ．初期値を摂動させたい場合，RandomReal[]を利用するとよい．
もちろん，NSolveやFindRootは連立方程式も解くことができる．
FindRoot を連立方程式にも用いる場合，Mathematica がちゃんと解釈できるように 初期値を与える必要があることに注意．たとえばf(x) = (f1(x), f2(x), f3(x))という関 数列とx= (x1, x2, x3)というベクトルがありf(x) = 0 という3本の方程式を解きたい とする．初期値をx1 = x2 = x3 = 0としたい場合は，以下のように3 × 2行列の初期値 を書く:
FindRoot[f[x] == 0, Transpose@{x,ConstantArray[0,Length@x]}]
初期値のほうは単に{x,0}と書くと {{x1, x2, x3}, 0} と解釈されてしまうので注意． ConstantArrayの部分はTableやAppendなどを使えばより複雑な初期値を設定するこ とができる．
3.2.4 FindMinimum
極値を求めるには FindMinimumを使うとよい．FindRootと違って，制約条件を明示 的に記述できて便利だ．ただ，FindRootよりは遅い．制約条件や探索範囲の部分を行列 や関数を用いて複雑化したり他の関数とネストしたりするとうまく動かないこともある． 適当な方程式f (x) = g(x)の解xを求めるときに，FindRoot[f[x] ==g[x],{x,0}] などとする代わりに，FindMinimum[(f[x]-g[x])^2, x]として最小化問題に変換する こともできる．連立方程式に拡張することもでき，f (x) = g(x) = 0 を求めたい場合， minxf (x)2+ g(x)2 とすると，最小化問題の解は連立方程式の解になる．FindMinimum だと制約付きのプログラムにも容易に拡張できるので，FindRootよりも幅が広がる．た だし実践上は，f やgの形状や定義域に注意も必要である．FindRootよりも誤差が大き くなるケースがある気もする．たとえばlocal minimaに収束してしまうと，f (x) = g(x) となる解には至らないことが考えられる．また，FindRootよりも計算が遅くなりがちで ある．停止条件をゆるめに設定したFindMinimumで解を求めたら，それを初期値として 改めてFindRootを回すといい気がする．
MaximizeやMinimizeと，その数値版NMaximize, NMinimizeもある．
3.2.5 LinearSolve
行列方程式，たとえばAx = bの解x = A−1bを求めるときにはLinearSolveを使う． 御利益がいかほどなのかは素人にはなかなか想像がつかないが，数値計算の世界では逆行
列を直接計算しないのが常識である．計算の誤差はゴキブリみたいに増えていくものだ． 固有値はEigenvaluesで，固有ベクトルはEigenvectorsで求めることができる．
3.2.6 DSolve, RSolve
微分方程式はDSolveで，差分方程式はRSolveで解く．たとえば，CES関数における 代替の弾力性が
σ = −f
′(x)[f (x) − xf′(x)] xf′′(x)f (x) , という微分方程式の解に対応していることを確認してみよう．
f[x] /. DSolve[sigma == -f’[x] (f[x] - x f’[x])/(x f’’[x] f[x]), f[x], x] // Simplify;
Simplify[% /. x -> (X^(sigma/(-1 + sigma)))^(1 - sigma), X > 0 && sigma > 0 && sigma < 1] ;
Simplify[% /. X -> (x^(1/(1 - sigma)))^((-1 + sigma)/sigma), x > 0 && sigma > 0 && sigma < 1]
TraditionalForm@%
ここで二回出てくるSimplifyでは適当に変数変換して簡略化を手助けしている．結果的 にはf (x)は次のように整理される:
{c2(c1 + xσ−1σ )σσ−1}
数値的に微分方程式を解くのはNDSolveで，第4.3章で実際の用例を簡単に紹介して いる．
3.3 図を書く
完成したグラフは，EPSやPDFファイルとして保存しておくとLATEXに持って行き やすい．該当のグラフを右クリック→形式を選択して画像を保存(A)，や，Exportコマ ンドなどを使えばよい(第3.3.6章)．
3.3.1 Plot, Plot3D, ListPlot, Show
作図コマンドPlotは直感的にわかりやすく，解析をよく助けてくれる．PlotLegends で凡例をつけたり，EpilogやPrologで低水準な図を重ねたり，という小技を覚えてお
くとプレゼンに使いやすい．
3次元のグラフを作るPlot3Dは，極値の大域的な特徴を目視確認することができ，ご つい分析をしている気分に浸ることもできるが，さほど生産的ではなくあまり使わない． 3D円グラフを使わないのと似たり寄ったりだ．
Plotで描写するオブジェクトが複雑な場合，うまく機能しないことがある．こういう 場合は，ListPlotで代替すると早い．
複数のグラフを重ねたり並べたりしたい場合は，ShowやGridやGraphicsGridを使 う．しかし，実践上はnotebook上できれいに描画したところで自己満足に終わることが 多く，役に立つ場面は少ない．
フォントを変えたい場合は，BaseStyle -> {FontFamily -> "Helvetica", 14}の ようにオプションをつければよい．
3.3.2 Manipulate
パラメータを動かすとグラフがどう変化するか，といった比較静学のようなことをイン タラクティブにできるのがManipulateである．たとえば
Manipulate[
Plot[{a x^2 + x + 1, - 2 a x^2 + x}, {x, - 5, 5}], {{a, 0}, -3, 3}]
のように，Plotと組み合わせてバーをぐりぐりいじる．発見的に解の性質を考えるのに 便利である．PlotだけでなくContourPlotなどと組み合わせて使うことも多い．CDF エキスポートすれば，Mathematica を持っていない人にこの面白い体験をお裾分けで きる．
問題によっては初期値を適切に設定しないとエラーで止まる．
3.3.3 ContourPlot, DensityPlot
解が陰関数でしか書けないことは珍しくない．こういうときでも解を絵的に確認したい 場合，等高線プロットContourPlotを使う．たとえば，適当に作図すると次の通り: ContourPlot[1/3 == x^2 + ( 5 y/4 - Sqrt[Abs[x]])^2,
{x, -1, 1}, {y, -1, 1}, ContourStyle -> Pink] ContourPlot[{x^2 + y^2 == 1,
x^(1/2) y^(1/2) == 1/Sqrt[2],
y == - x + Sqrt@2},
{x, 0, 1.2}, {y, 0, 1.2}, Epilog -> {Point[{1/Sqrt@2, 1/Sqrt@2}]}] 以下の図のような出力が得られる．
-1.0 -0.5 0.0 0.5 1.0
-1.0 -0.5 0.0 0.5 1.0
0.0 0.2 0.4 0.6 0.8 1.0 1.2
0.0 0.2 0.4 0.6 0.8 1.0 1.2
ContourPlotと似たものに，密度プロット DensityPlotがある．個人的には使った ことがないが，数学の達人たちの論文でたまに見る．
3.3.4 RegionPlot
ある不等式を満たすパラメータ(a, b)の組が描く領域を知りたいとする．ContourPlot
で(a, b)平面上に方程式を書いてもよいが，こういう場合はRegionPlotのほうが便利
だ．たとえば，
regPlot1 = RegionPlot[{a^2 + b^2 <= 1}, {a, 0, 1.2}, {b, 0, 1.2}]; regPlot2 = RegionPlot[{a^(1/2) b^(1/2) >= 1/Sqrt[2]},
{a, 0, 1.2}, {b, 0, 1.2}, PlotStyle -> ColorData[97, "ColorList"][[2]]]; Show[regPlot1, regPlot2]
のようにする(下図)．ちなみにColorData[97, "ColorList"]では，Mathematica 10 デフォルトのカラーテーマを指定している．
なお，RegionPlotはこの例のように2つに分けて書く必要は必ずしもなく，&&(and) などの論理演算子を使えばもう少しシンプルになる．実践上も，複数の不等式を同時に満 たすパラメータの範囲を知りたいことも少なくない．RegionPlotはそれなりに使用頻度
0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0
0.2 0.4 0.6 0.8 1.0 1.2
が高く，頭の片隅においておくべき関数の一つである．
ContourPlotやRegionPlotは，グラフが非連続になるようなポイントや解が一意に 得られないとき，描画が思うように行かないことがある．こうした場合，MaxRecursion
やPlotPointsオプションをいじると，よりきれいな図が得られることがある．(ただし
ファイルサイズが少し重くなる．) ただし，MaxRecursionやPlotPointsをいじるより も，プロットすべき関数を書き直すほうがよいことも多い．
version 10 で ，ImplicitRegion や ParametricRegion な ど を 使 っ て 複 雑 な 領 域 を 定 義 す る こ と が で き る よ う に な っ た ．領 域 の 和 や 積 や 差 は RegionUnion, RegionIntersection, RegionDifferenceを使えばよいし，より論理演算をがっつり 使いたい場合はBooleanRegionが使える．
3.3.5 ParametricPlot
x座標もy座標も何らかの関数になるようなグラフをプロットしたいとする．たとえば (x, y) = (f (z), g(z))という座標を持つグラフをz を動かしながら(x, y)平面に図示した いとする．こういう場合はParametricPlotを使う．
凸包(convex hull)を図示したい場合，ParametricPlotを使ってもよいが，より直接
的にConvexHullMeshを用いるとよい．たとえばくり返しゲームのフォーク定理を説明
するときに出てくるようなfeasible setを書きたいとき， payoff = {{0, 0}, {2, 2}, {3, -1}, {-1, 3}}; ConvexHullMesh[payoff, Axes -> True]
と利得行列を直接放り込めばよい(下図)．
-1 1 2 3
-1 1 2 3
領域の境界部分に関心がある場合，RegionBoundaryを用いればよい．
3.3.6 プレゼンテーション用画像の作成
私がプレゼンにMathematicaの図を使う場合にどうやっているか紹介する．
私の場合，講義などのカジュアルなプレゼンや自分用メモには PowerPointを，意識 高そうなオーディエンス向けプレゼンにはKeynoteを，学術用のプレゼンでは Beamer を使っている．画像の種類はpdfかepsかpngである．TEXで作る文書・レジュメには pdf(またはeps)を用いるので，たいていの場合pdfにしている．(tikzで扱えるとよいの だがやり方がよくわからない．)
Mathematicaからpngなど画像を出力するには， 1. Exportコマンドを使う
2. ノートブック右端の]を右クリック → 「選択範囲の形式保存…(L)」
3. ノートブック右端の]を右クリック → 「選択範囲の印刷…(R)」→ PDFにする 4. 画像を右クリック→ 「形式を選択してグラフィックスを保存」
5. 画像を右クリック → 「グラフィックスをコピー (C)」 → Adobe Illustrator に ペースト → PowerPointに図として貼り付け
といった手段を用いる．(フォントの埋め込みなど細かい仕様はよく知らない．) 5番目の 右クリック→コピーは直感的で楽である．凡例もコピーしたい場合はマウスをドラッグし 凡例ごと選択してから右クリックする．
epsやpdfにした場合，背景など余計な要素がくっついてくることもある．こういう場 合はAdobe Illustratorを用いて，グループ解除→クリッピングマスクを解除→余計な 要素を削除，と加工することもある．(もちろん結果を恣意的に操作する意図はなく，単 に見栄えを整えるためである．悪用は厳禁である．) Illustratorではさらにフォントやカ ラーを調整することも多い．
Manipulateなどのインタラクティブな画像を出力して使いたい場合は，「ファイル(F)
→ CDFエキスポート(D)」を使う．ウェブサイトなどにあらかじめアップして，アクセ
スすれば使えるようにしておくと便利だ．私の経験上，インタラクティブなプレゼンは， 非専門家とコミュニケートするときにとりわけ効果的となる．が使う機会は少ない．
3.4 その他
Mathematicaは「リスト」という概念で，何かをまとめて扱う．ベクトルや行列もリス
トである．リストは入れ子にもできる．リストの扱いに慣れることが上達の鍵だと思う．
PythonやRと併用すると混乱するポイントでもある．
たいていの関数は，リストの個別の要素に適用される．たとえば Log[{a,b,c}] は， {Log[a], Log[b], Log[c]}を返す．ベクトルでなく行列でも同様である．
3.4.1 Clear
定義した変数や関数を一度リセットしたいときは，Clearを使う．notebookの一番最 初に実行するセルの最初の行に
Clear["Global‘*"]
と入れておくと，まっさらにして一からスタートできる．
3.4.2 Table
分野問わず Tableは重要な関数であり，Mathematicaを使いこなす上での鍵となる． 何らかの規則性があるベクトル・行列・テンソルを生成するのに使う．MatlabやRのよ
うなMathematica以外のソフトを学ぶと，Tableの活用法が見えてきやすい．実際の用
例は第5章でも取り上げる．
簡単な例として，この原稿を執筆しているのは 2016年であるが，2016 = 33 + 43 + 53+ 63+ 73+ 83+ 93,である．Tableを用いてこれを確認するには，
Total@Table[i^3, {i, 3, 9}]
とすればよい．(この例のような単純な数列の和ならSum[i^3, {i, 3, 9}]のほうがよ りよい．)
Tableと類似の関数に，RangeやArrayがある．RangeはTableの簡易版で，ちょっ とした等差数列を作るだけならRangeで十分だ．
Tableは評価のタイミングが少し複雑で，他の関数と組み合わせるときにうまくいかな
いときもある．(Table[expr, spec]は, たとえばspecの部分に局所化できない何かを 入れているとうまくいかない．)こうした場合，Tableの代わりに純関数を用いて書くと うまくいくことがある．(詳しい仕組みはよくわからない．)
3.4.3 Transpose, Tr
行列を操作するときに転置 Transpose は頻出である．「Esc tr Esc」でもよいが，
Transposeのほうが読みやすいようにも思う．好みで．
トレースはTrになる．似ているので注意．
3.4.4 Integer, Rational, Real, Complex
数のクラスは Integer, Rational, Real, Complexの4つである．意味は読んで字の ごとし．
虚数の実部・虚部を取り出すには，それぞれRe, Imを使う．もちろん普通の実数なら Im[x] = 0となる．
3.4.5 AllTrue, AnyTrue, NoneTrue, Cases, DeleteCases, パターン
生成したリスト(主にベクトルや行列)に，虚数や負の数が含まれていないか知りたい ことは少なくない．このように，何らかのパターンにマッチする要素を操作したいときは いくつかやり方がある．version 10から実装されたAllTrue, AnyTrue, NoneTrue, が直 感的には使いやすいと思う．
たとえばベクトル xの要素に厳密に正の値が含まれていないこと(x ≪ 0)を確認した い場合，
AllTrue[x, NonPositive] NoneTrue[x, Positive]
などでTrue or Falseを見ればよい．
その他，何らかのパターンにマッチする要素を操作したいときはCasesやDeleteCases を使うこともできる．欠損値を省きたい(R でいうna.omit())ときは，Delete も利用 できる．
虚数を取り除きたい場合，たとえば {1, 3, 4 + I, -1.5};
DeleteCases[%, _Complex]
とすると，虚数を含む第3要素だけが消え去ったベクトルが生成される． 負の値を消したい場合，Signを活用してもよいし，Ifを使って
Table[RandomReal[] - 1/2, {5}] (* 判別したいリスト *)
If[# >= 0, #, negative] & /@ % (* もし負だったらnegativeと出力 *) DeleteCases[%, negative]
のようにしてもよい．リストに虚数が混じっているとIfのところでエラーが生じるの で，その可能性はあらかじめつぶしておく．
消さずにその位置だけ知りたい場合，Positionを使えば良い．リストの次元を変えず におけるので便利だ．
リストxが空である(x = {})かどうか知りたいときは MatchQ[x, {}]
とすればTrue or Falseで判定してくれる．またxが空っぽであればLength@xとしたと き0と返ってくる．
_Integer は整数かどうかを判別するのに使うパターンだが，0 は整数で0.は整数で
はないと認識されるので注意(第2.1.9章)．
CasesやIfなどは，複雑な関数を組むときに活躍するので頭の片隅においておきたい．
3.4.6 反復，選択
他のプログラム言語でも使えるようなループ関数や条件関数，For, While, Do, If, Which, Switch, はMathematicaでも使える．使い方は普通のプログラミングのテキス トを読むか，適当な例を探すべし．仕組みを押さえればMathematicaの運用能力も向上 する．なお，Tableなどの組み込み関数で代替できる場合はそのほうが早いと思う．
使えると言っても，Mathematicaはこの手の関数をごりごり使う用途のソフトではな いと思う．反復するだけなら他のソフトでもできるしその方がシンプルで高速になる気が
する(Pythonなど)．
3.4.7 Module
自作で関数を作る際，その手続きの中だけで使いたい変数 (局所変数) があるときは Moduleを用いる．関数を f[x_]:=exprなどと定義するよりも応用範囲が広い．プログ ラミングの素人(含む私)はつい避けてしまいがちだが，手続きが長く複雑な数値計算を するときなどには活用したほうがよいと思う．
類似の関数にBlockやWithがある．私はあまり使ったことがない．
ModuleとBlockとの違いは，Blockがシンボルに割り当てられた「値」を局所化する のに対して，Moduleは局所変数の「名前」を局所的に扱う，という点である．Moduleと Withとの違いは，Withは局所「定数」を使用するのに対して，Moduleで使う局所「変 数」を使うことである．前者は値が1回評価されたらそのまま動かないが，後者は動きえ るのである．
3.4.8 FixedPoint, Nest, NestList
不動点を見つけたい場合やその他適当な反復(iteration)アルゴリズムを収束するまで 走らせたい場合，FixedPointが使える(第4.5.2章)．途中経過の履歴を残したいときは
FixedPointListを使う．しかし，工夫しないと計算時間がかかるケースも多い．最大ス
テップ数や停止幅を緩やかに定めておくとよい．
単に適当な回数同じ関数を適用するときにはNest，またその結果を毎回を格納してい くリストを作りたい場合はNestListを使う．要は，(f ◦ f ◦ · · · ◦ f)(x)のようなn回分 の合成関数を計算するものである．不動点問題以外でも，もろもろのシミュレーションで 活躍する．
3.4.9 計算時間が気になる
計算時間がどのぐらいなのか気になるときは，測れば良い．Timingを使ってもよいが， 長いプログラムでは次のように書けば何秒かかったか教えてくれる:
timebegin = SessionTime[]; (…計算時間を知りたいプログラム…)
"time used"
SessionTime[] - timebegin
しばしば，計算時間が長引くのは Mathematicaやマシンのせいではなく，プログラ
マー側の責任である．
3.4.10 音を鳴らす
評価が終わったタイミングが気になるときは，音が鳴るように設定しておけばよい．た とえばラ(A)の音(440Hz)を鳴らすには，
EmitSound@Sound@SoundNote[9]
でよい．9を省略するとド(C)の音になる．Soundのオプションをいじればピアノ以外の 音色も選べる．正弦波が好きなアヴァンギャルドな方は，
EmitSound@Play[Sin[2 Pi 440 t], {t, 0, 0.5}]
として0.5秒間440Hzを堪能してください．(個人的にはもっと高い周波数が好き)
ビープ音では味気ない，という人は好きな音楽ファイル (MIDIや WAV) を用意し， Import["E:\\Dropbox\\hogehoge.mid"]などと読み込み，計算の待ち時間中に再生ボ タンを押せば良い．こうすればマリオやドラクエの効果音を計算完了時に流すことがで きる．
3.4.11 組込み関数の中身が知りたい…が
Mathematicaが標準装備している関数の具体的な中身は非公開であり，リバース・エ
ンジニアリングは難しい．
4 実践例編
実践例を教育的に紹介する．ミクロなら第4.1, 4.2章，マクロなら第4.3章，あたりを 参考に．
4.1 単純な例 : Stackelberg model
notebookは，共著者が見ても何を計算しているのかがわかるように書くとよい．
たとえば，線形需要・限界費用一定の素朴なStackelberg競争を解く場合，私なら次の ように書く．
(* defining payoffs *) p = a - b (x1 + x2);
profit = {(p - c) x1, (p - c) x2};
"stage 2: firm 1 chooses x1 given x2"
foc1 = Solve[D[profit[[1]], x1] == 0, x1] // Simplify // Flatten
"stage 1: firm 2 chooses x2 with x1(x2)"
foc2 = Solve[D[profit[[2]] /. foc1, x2] == 0, x2] // Simplify // Flatten
"equilibrium outputs"
focs = Flatten@{foc1 /. foc2, foc2} // Simplify
"equilibrium price" p /. focs // Simplify
"equilibrium payoffs"
eqprofit = profit /. focs // Simplify
ただし，二階条件のチェックは省略している．場合分けをすっ飛ばして内点解だけに注 目してさっさと解きたい場合は，シュタッケルベルグに限らずだいたいこのようなコード になる．Dで微分し，Solveで解いて，Simplifyで整理し，Flattenで余計な括弧を外 す，/.で得た式を代入する，という流れは汎用性が高い．
その他，解く手順をスムーズにする，解いた結果を適当な変数に格納する，変数名はそ れなりに意味が通るようにしつつ，ラフな英語で適宜補注を入れる，必要に応じてベク トルを使い整理する，といったところが主なポイントとなる．前置修飾・後置修飾の順番 は，結果に影響がなければ気分次第で．
こうして得たfocsを評価すると，均衡生産量が {
x1 → a − c4b , x2 → a − c2b }
と出力され，eqprofitを評価すると均衡利潤が以下のように得られる: { (a − c)2
16b ,
(a − c)2 8b
}
適当なパラメーターを与えて，均衡を図示してみよう． para = {a -> 1, c -> 1/10, b -> 1};
bresponse = x1 /. foc1 /. para; focsPara = {x1, x2} /. focs /. para;
profit2 = {profit[[2]], eqprofit[[2]]} /. para;
ContourPlot[{bresponse == x1, profit2[[1]] == profit2[[2]], profit2[[1]] == 0.5 profit2[[2]],
profit2[[1]] == 1.5 profit2[[2]]}, {x1, 0, 1}, {x2, 0, 1}, Epilog -> {PointSize[0.02], Point[focsPara]},
FrameLabel -> {"x1", "x2"}];
Show[%, Graphics[{Text[ Style["Equilibrium", 14], {0.65, 0.22}], Arrow[{{0.6, 0.25}, {0.25, 0.43}}]
}]]
Equilibrium
0.0 0.2 0.4 0.6 0.8 1.0
0.0 0.2 0.4 0.6 0.8 1.0
x1
x2
後手の最適反応曲線と，先手の等利潤曲線と，均衡点がわかるようになっている．おま け的な部分は数値を手動で入力しているのであまり美しくないが，およそこのような要領 でグラフを補助的に作成していく．前頁のような出力を得る
4.2 変数変換の例 : Cournot duopoly
この章では，limited domain problemへの対処について触れる．そしてJudd (1998) のchapter 5.4にある例をMathematicaで再現する.*1
limited domain problemは文字通り，変数の定義域が限定的な状況でどう数値計算す
るか，を問題にする．経済学で扱う変数の多くは非負であり，解きたい問題 f (x) = 0 がある状況で，x < 0のとき f がうまく定義できないことがある．たとえば効用関数 u(x) = log(x)やu(x) = √x = x1/2などだ．一方で，特に何も工夫せずニュートン法を 回すとR全域を探索してしまうので，しばしf (x)が虚数になってしまい失敗する．
MathematicaでFindRoot[f[x]==0, {x, x0, xmin, xmax}]と書くとxの定義域 を[xmin, xmax]に限定することができる．しかしこの場合，定義域の外側に飛び出た瞬間 に試行を終了してしまい，f (x) = 0という解にたどり着くことがほぼ期待できない．
いくつか対処法はありえる．ここでは定義域をRになるように変数変換をしてみよう． クールノー競争の均衡を求めたい．2財(Y とZ)あり，限界費用をそれぞれcY, cZ と おく．消費者の準線形効用関数を以下のようにおく(M はニュメレール):
U (Y, Z) = u(Y, Z) + M = (1 + Yα + Zα)η/α + M.
Y財を作る企業の利潤はΠY(Y, Z) = (u1− cY)Y , Z財を作る企業の利潤はΠZ(Y, Z) = (u1− cZ)Z である．u1, u2 はそれぞれ偏微分を表す．求めたい均衡は次の連立方程式の 解である:
ΠY1(Y, Z) = ΠZ2(Y, Z) = 0.
Limited domain problemの問題は，Y, Z ∈ R+であることだ．そこで， y = log Y, z = log Z,
と定義した変数を用いて，
ΠY1 (ey, ez) = ΠZ2(ey, ez) = 0,
と解きたい問題を書き直す．y, z ∈ Rについてニュートン法を走らせるのである．
*1Judd, Kenneth L. Numerical methods in economics. MIT press, 1998.
パラメーターをα = 0.999, η = 0.2, cY = 0.07, cZ = 0.08とする．初期値 (y0, z0) = (−2, −2) から探索すると答えは (y∗, z∗) = (−0.137, −0.576), あるいは (Y∗, Z∗) = (0.87, 0.56)に至る．コードは以下の通り
(* define parameters and objective functions *)
para = {alpha -> 0.999, eta -> 0.2 , cY -> 0.07, cZ -> 0.08} U[Y_, Z_] := (1 + Y^alpha + Z^alpha)^(eta/alpha) + numeraire profitY[Y_, Z_] := Y (D[U[Y, Z], Y] - cY) /. para
profitZ[Y_, Z_] := Z (D[U[Y, Z], Z] - cZ) /. para
(* reaction functions and transformation *)
reactionY = D[profitY[Y, Z], Y] /. {Y -> Exp@y, Z -> Exp@z} // Simplify reactionZ =
D[profitZ[Y, Z], Z] /. {Y -> Exp@y, Z -> Exp@z} // Simplify
(* solutions *) y0 = -2; z0 = -2;
solution = FindRoot[{reactionY == 0, reactionZ == 0}, {y, y0}, {z, z0}] Exp@{y, z} /. solution;
solution2 = {Y -> %[[1]], Z -> %[[2]]}
解に収束するプロセス(5ステップかかる)と合わせて最適反応関数を図示した: (* monitoring *)
{result, steps} = Reap[FindRoot[{reactionY == 0, reactionZ == 0}, {y, y0}, {z, z0}, StepMonitor :> Sow[{y, z}]]]
(* reaction curves. See also Judd, Fig 5.7 *) ContourPlot[{reactionY == 0, reactionZ == 0},
{y, -2, 0.5}, {z, -2, 0.5}, FrameLabel -> {y, z},
Epilog -> {Red, PointSize[Medium], Map[Point, Prepend[steps[[1]], {y0, z0}] ]} ]
-2.0 -1.5 -1.0 -0.5 0.0 0.5 -2.0
-1.5 -1.0 -0.5 0.0 0.5
y
z
4.3 微分方程式 : Solow model
学部レベルのソロー・モデルを2通りのやり方で動かしてみる．一つは線形近似，もう 一つはNDSolveを使った数値計算，である．
一人あたり資本ストックをkt, 一人あたり産出量をyt = Aktα, とする．全要素生産性 A，人口成長率n，貯蓄率s，資本減耗率δ, 資本分配率α,初期の資本ストックk0, をす べて外生一定とする．モデルはおなじみの非線形微分方程式に帰着する:
˙kt = g(kt) ≡ sAkαt − (δ + n)kt, (1)
定常状態における資本ストックをk∗ とおくと，次のように整理できる:
k∗ =
( sA δ + n
)1−α1 . (1)式を定常状態の周りで線形近似すると，
˙kt ≃ g′(k∗)(kt − k∗),
である．˙kt = kt− k∗ を使うと簡単な自励系の表現となり，解はkt ≃ k∗+ Ceg′(k∗)t,と書 ける．ただしCは初期条件から決まる任意の定数であり，t = 0を代入すればC = k0−k∗
と解ける．今g′(k∗) = −(1 − α)(δ + n)であるから，
kt ≃ k∗+ (k0− k∗)e−(1−α)(δ+n)t.
ここまでの計算をMathematicaで簡単に確かめよう．
(* right hand side of the law of motion *) g[k_] := s A k^alpha - (delta + n) k
(* k at steady state *)
kStar = Solve[g[k]== 0, k] // Simplify // Flatten (* growth rate of (k - kStar) *)
growth1 =
Simplify[D[g[k], k] /. kStar ,
1 > alpha > 0 && A > 0 && s > 0 && n + delta > 0] これを実行すると次のような出力を得る:
{ k →
( As
delta + n
)1−alpha1 }
(alpha − 1)(delta + n)
得られた資本ストックを用いて，生産性の恒常ショックが与える影響を図にして確かめ てみよう．0期にA = 1 であったものが，t = 1, 2, . . . , 24までA = 1.5に上昇したとす る．残りのパラメーターは次のように適当に設定する:
para1 = {alpha -> 1/3, delta -> 1/20, n -> 1/20 , s -> 4/10, A -> 1}; para2 = {alpha -> 1/3, delta -> 1/20, n -> 1/20 , s -> 4/10, A -> 1.5};
初期のパラメーター(para1)下で，定常状態における資本ストックを計算するとk∗ = 8 となる．初期の資本ストックをk0 = k∗ = 8としよう(k0Star)．1期目以降はパラメー ターに変化があり(para2)，資本ストックはkLA[t]という関数に従って動くことになる．
(* initial value, at steady state *) k0Star = k /. kStar /. para1
(* derived path of capital: a linear approximation *)
kLA[t_] := k + (k0 - k) Exp[growth1 t] /. kStar /. para2 /. k0 -> k0Star Plot[kLA[t], {t, 0, 24}, AxesLabel -> {t, k[t]}]
(*---illustration---*)
Plot[If[t < 0, k0Star, kLA[t]], {t, -6, 24}, AxesLabel -> {t, k[t]}]
最後の行では簡単に作図して確かめている(次図)．もちろん産出量や消費・貯蓄など他 の内生変数も同様に作図することができる．
-5 5 10 15 20
t 9
10 11 12 13 k(t)
次に，NDSolveを使ってMathematicaに一連の微分方程式を数値計算してもらおう． kmotion = k’[t] - g[k[t]] /. para2
kNDSolve = NDSolve[{kmotion == 0, k[0] == k0Star}, {k[t]}, {t, 0, 24}] (* evaluate the interpolating function to visualize *)
kNDSolve2[t_] := Evaluate[kNDSolve[[1, 1, 2]]] (* ---illustration--- *)
Plot[{kNDSolve2[t], kLA[t]}, {t, 0, 24}, PlotLegends -> {"NDSolve", "Linearize"}, AxesLabel -> {t, k[t]}]
グラフを比較する(下図)と，線形近似で解いたパスと，数値計算によるパスは若干の 誤差があることがわかる．線形近似のほうがexactな解との乖離は大きくなるだろう．
5 10 15 20
t 9
10 11 12 13 k(t)
NDSolve Linearize
最後に，(1)式を使い，初期値 k0 を与えてそこから定常状態k∗ までのkt の移行過程 を追いかけてみよう．(1)式の離散時間版，kt+1 = g(kt) + kt, を念頭に置き，k0 = 4か らスタートした場合とk0 = 6からスタートした場合について，t = 10までの資本ストッ
クの推移を図示すると次のようになる: tT = 10;
x = Table[k, {k, 0,tT}];
kevolve = NestList[Evaluate[g[#] + # /.para1] &, k, tT] ; {Transpose@{x, Evaluate[kevolve /. k -> 4]},
Transpose@{x, Evaluate[kevolve /. k -> 6]}};
ListPlot[%, Joined -> True, PlotLegends -> {"k_0=4", "k_0=6"}, AxesLabel -> {"t", "k_t"}]
こうしてabsolute convergenceのプロセスが視覚的に確認できる．なおg(k) + k をく り返し計算している3行目(kevolve)では，Evaluateで数値化している．もしこれが なければ，毎回gの右辺を呼び出して1/3乗を含む多項式の1/3乗を含む多項式の1/3 乗…などを無限精度で求めてしまい，計算時間がとてつもないことになるようだ．そうで なくてもtが大きい場合は遅いので注意．
4.4 agent-based model の例 : Schelling’s segregation model
PythonやJuliaでおなじみ，Sargent and StachurskiのQuantEcon掲載の例をMath- ematicaで再現してみる.*2 T.C. Schelling (1969 AER)の人種による分離居住パターン の発生をシミュレートしてみたい．
詳しい解説はQuantEconを参照してもらうとして，主な設定は次の通りだ． プレイヤー 2タイプの人がいる．ここでは250人ずつ存在．
立地 居住する場所を2次元平面(x, y) ∈ [0, 1]2から選ぶ．
戦略 プレイヤーは，近隣に同じタイプの者がどれだけいるか，を気にしている．特に， 最も近い10人のうち，5人以上が自分と同じタイプなら「happy」とする．
• happyなら，立地はそのまま．
• happyでなければ，別のhappyになれそうな場所に移動する．
こうした簡単なルールから，マクロレベルで見ると秩序あるパターンが生み出される． さっそくコードを書いて確認しよう．まずはパラメーターの設定と，各人のタイプを表す ベクトル(type)，初期の立地(座標locationX, locationY)を設定する．
*2Thomas J. Sargent and John Stachurski (2017) QuantEcon.lectures-python3 PDF.
