Rのグラフに回帰分析の結果である数式を書き込む

データをプロットしたグラフに、推定された予測値や、推定されたパラメーターを数式ごと書き込みたいときはありませんか？
そう難しい話では無いのですが、幾つかコツがあるので例を作ってみました。

1. データセット

グラフを描くにはデータセットが要るので、今回は以下の数字を

# 奈良女子大学のページにあった実験観測データを拝借
t <- c(0,5,10,15,20,25,30) # 時間
Ts <- c(56.5, 32.1, 21.9, 17.0, 14.6, 11.9, 10.9) # 物体の温度
Tm <- 0 # 外気温（冷却材は氷でゼロ度だが、室温はもっと高い模様。なお、9.660583度を仮定すると、もっとも当てはまりが良い）

2. 推定を行う

推定を行い、予測値を計算してみましょう。モデルの制約によりpredict関数は使えません。

# ニュートンの冷却の法則を想定し、定数を推定する
dT <- (Ts - Tm)/(Ts[1]- Tm)
# OLSで推定。
r <- lm(log(dT) ~ t + 0)
# 予測値を計算
prdct <- Tm + (Ts[1] - Tm)*exp(coef(r)[1]*t)

3. プロットする座標の範囲を決める

グラフを重ねる前に縦横の座標の範囲を決めておかないと、ズレます*1。

xlim <- c(0, max(t))
ylim <- c(0, 100)

4. データをプロットする

まずは観測値をプロットしてみます。

par(mar=c(5, 3.5, 1, 1), mgp=c(2, 1, 0))
plot(Ts ~ t, xlab="経過時間", ylab="", axes=FALSE, xlim=xlim, ylim=ylim)
x_at <- seq(xlim[1], xlim[2], 10)
y_at <- seq(ylim[1], ylim[2], 25)
x_labels <- sprintf("%d分", x_at)
y_labels <- sprintf("%d℃", y_at)
axis(1, at=x_at, labels=x_labels, pos=ylim[1])
axis(2, at=y_at, labels=y_labels, pos=xlim[1], las=2)

5. 推定値を重ね書きする

推定値をグラフに書き込みます。ここはRの紹介で良く触れられますね。

par(new=T)
plot(prdct ~ t, type="l", xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim, lty=2)

なお、こんな手間隙かけずにlines(t, prdct, lty=2)するべきだろと突っ込まれました。おっしゃる通り。
直線なのが気に入らない人は、plotの代わりに以下のようにcurveを使ってください。

par(new=T)
curve(Tm + (Ts[1] - Tm)*exp(coef(r)[1]*t), xname="t", xlab="", ylab="", axes=FALSE, xlim=xlim, ylim=ylim, lty=2)

6. 推定された式を凡例として表示する

まず、数式をsprintf関数で作って、parse関数で表現式に変えます。expression関数は、代入が手間なので使いません。使える数式の記号はhelp plotmathをすると一覧が見られます。

expr <- parse(text = sprintf("%.3f + (%.3f-%.3f)*e^{%1.3f*t}", Tm, Ts[1], Tm, coef(r)[1]))

次に、legend関数で凡例を表示します。推定値はlty、観測値はpchになりますが、使わない方には-1を入れていることに注意してください。

legend("topright", lty=c(2, -1), pch=c(-1, 1), legend=c(expr, "観測値"), bty="n")

7. データの出所を明記する

データを拝借しているので、出所を明記しておきましょう。

mtext("出所) 奈良女子大 - 「第11回の授業　(07/07/09)」の実験データから編集" , side=1, line=3.5, adj=0)

8. 完成図

相関係数などを書いても良かったかも知れません。

9. 補足

文字が化ける場合は、以下のようにpar()関数で指定するグラフィックス・パラメーターでフォントを指定してください。フォントはMeiryoを指定していますが、お使いのOSにあわせましょう。

windowsFonts(GraphFont = windowsFont("Meiryo"))
par(mar=c(5, 3.5, 1, 1), mgp=c(2, 1, 0), family="GraphFont")

画面のグラフを保存したい場合は、以下のような感じで。

dev.copy(png, "ニュートンの冷却の法則の実験結果.png", width=640, height=480, bg="white")
dev.off() # 画面のグラフを消す
dev.off() # PNGファイルを閉じる

*1:推定値のプロットにpar(new=T)とplot()を使うからズレるのであって、低水準関数のlines()を使えば座標の範囲を決めなくてもいいと言うツッコミがありましたが、ここを決めておかないと軸目盛などの位置もズレると言うか、あわせるのが困難になるので、やはり決めておいた方がいいです。

Rのグラフ中文字列の上下左右の寄せ

text関数でプロット領域に文字を書き込むことができるわけですが、指定座標と文字位置の関係を示すadjパラメーターの意味を整理してみます。
グラフの中心に表示座標をとり、adjだけ色々と変えてみたグラフが以下です。

• c(1, 0)だと、文字列の右下に表示座標が来ます。
• c(0, 0)だと、文字列の左下に表示座標が来ます。
• c(0, 1)だと、文字列の左上に表示座標が来ます。
• c(1, 1)だと、文字列の右上に表示座標が来ます。

0.5を入れておくと中央になるわけですね。1を超えた値や、マイナスの値も有効です。
mtext関数のadjはベクトルではなく、値一つのスカラー値になるので注意してください。上下寄せを指定できません。

Rのグラフで余白の調整

Rのプロットはデフォルトでは余白が多いです。しかし、TeXの方でタイトルや注釈を入れる場合は、プロット領域を大きく取りたくなるので、この余白を削りたくなります。par関数のmar、mgpオプションで余白を調整できるので覚えておきましょう。

1. marオプション

marオプションは、下・左・上・右の余白を指定します。

デフォルトでは以下のようになっていますが、タイトルや注釈・出所の有無で調整した方が良いでしょう。

par(mar=c(5, 4, 4, 2) + 0.1)

単位はmexで行の高さです。

2. mgpオプション

mgpオプションは軸ラベル・軸メモリ・軸線の位置を指定します。

デフォルトでは以下のようになっていますが、軸ラベルの位置がプロット領域から遠すぎるので、c(2.5, 1, 0)ぐらいを指定した方が幸せになれると思います。

par(mgp=c(3, 1, 0))

単位はmexで行の高さです。

3. タイトル

plot関数でmainオプションを指定しないで、plot後にtitle関数でlineオプションをつけるのが良いようです。

title(main="グラフの表題", line=1)

Rで多体問題を数値解析

古典物理学と言えばニュートンの運動方程式で、何かの軌道をプロットするのが定番です。三体以上の問題は特殊なケースしか閉じた解が出てこないので、数値解析が行われています。銀河にある星の数をシミュレーションすると、スーパーコンピューターな数値演算能力が必要な世界になるわけですが、二個や三個であればPC上のRで計算できるので、試しにやってみましょう。

1. 計算する式

まずは物理の教科書に書いてありそうなニュートンの運動方程式の確認から*1。

が万有引力定数、が質量、が距離、は質点の数、は質点、が時点を表します。
力の成分をとに分解します。成分は省略。

両辺をで割って、、と置いて、常微分方程式に変形しましょう。

これをルンゲ＝クッタ法などで計算します。二体問題でも解析解を出すのには色々と苦労するので、手抜きな感じがたまりません。

2. ソースコード

N=3でも4でもいいのですが、N=2のソースコードを提示します。まずは演算部分。

# deSolveパッケージを使います
library("deSolve")


# 常微分方程式化したモデル
model <- function(t, X, params){
  with(as.list(params),{
    # 引数を分解
    x <- X[1:n]
    y <- X[(n+1):(2*n)]
    dx1 <- X[(2*n+1):(3*n)]
    dy1 <- X[(3*n+1):(4*n)]
    # 戻り値の領域を確保
    dx2 <- numeric(n)
    dy2 <- numeric(n)
    for(i in 1:n){
      r3 <- sqrt((x[-i]-x[i])^2+(y[-i]-y[i])^2)^3
      dx2[i] <- sum(-G*m[-i]*(x[-i]-x[i])/r3)
      dy2[i] <- sum(-G*m[-i]*(y[-i]-y[i])/r3)
    }
    list(c(dx1, dy1, dx2, dy2))
  })
}


# 計算のループ条件。0.1刻みで250まで。
times <- seq(0, 250, 0.1)
# 質点の数
n <- 2
# 質点の重さ
m <- c(500, 1)
# 初期座標と初期速度を以下の並び順で指定
# x[1], x[2], y[1], y[2], dx[1], dx[2], dy[1], dy[2]
X <- c(100, 400, 100, 400, 1, -1, 1, 0)
# 常微分方程式として計算する
out1 <- ode(X, times, model, parms = list(n=n, G=-9.80665, m=m), method="rk4")

演算結果は数値の羅列で分かりづらいので、グラフにしましょう。

# プロット範囲を指定
xlim <- c(0, 800)
ylim <- c(0, 800)


# N時点までのプロット
plot_as_of_N <- function(N){
  plot(out1[,2][1:N], out1[,2+n][1:N], xlim=xlim, ylim=ylim, type="l", xlab="x", ylab="y", lty=2)
  for(i in 2:n){
    par(new=T)
    plot(out1[,i+1][1:N], out1[,i+1+n][1:N], xlim=xlim, ylim=ylim, type="l", lty=i+1, xlab="", ylab="", main=paste("t=", round(N)))
  }
  for(i in 1:n){
    points(out1[N,i+1], out1[N,i+1+n], pch=i, cex=2)
  }
  legend("topleft", lty=seq(2, length(m)+1), legend=sprintf("m=%d", m), bty="n")
}


# 同時に2×2個をプロット
par(mfrow = c(2,2), oma = c(0, 0, 0, 0))
for(i in seq(1, length(out1[,1]), (length(out1[,1])-1)/3)){
  plot_as_of_N(i)
}

z軸を加えたり、質点の数を増やすのは難しく無いと思います。例えば以下のように三つのパラメーターをいじれば、三体問題化できます。

n <- 3
m <- c(500, 2, 1)
X <- c(100, 400, 800, 100, 400, 0, 1, -1, 0, 1, 1, 0)

衝突処理が無いので質点が近づきすぎると、どこかに飛んで行くのでパラメーターの設定には注意してください。

3. 計算結果

パラメーターに大きく依存するわけですが、上のコードだと下のようなグラフが描けます。最近は波ばかりプロットしていたので、グルグルしていて爽快ですね。

なおオイラー法などで計算すると、誤差が大きくなって全く違う図になります。

*1:詳しい説明は富久(2010)の第5章を参照

RでKdV方程式を数値解析でプロットしてみる

浅水波を表すKdV方程式は、今では可積分系として閉じた解が知られている非線形偏微分方程式ですが、これは意外に最近知られるようになった事で、二十世紀の中ごろにザブスキーとクルスカルと言う物理学者がシミュレーションをしてソリトンと言う孤立波を見つけ、そこから解析が大きく進むようになったそうです（武部(2008)）。当時を偲んでRで数値計算してみましょう。

1. KdV方程式の形状

まずは式の確認です。

が時間、が位置、が波の高さになります。
このKdV方程式は、例えば波動方程式と比べると厄介な形状をしています。三階微分が中にあり、数値演算のお決まりのパターンである常微分方程式にならないのです。

2. 数値演算が出来る形に変形する

どうしたら良いのかな・・・と思っていたら、Numerical Solution of the KdVと言うページに全てのノウハウが書いてありました。
まずは第二項を整理します。合成関数の微分を思い出すと何をしているのか分かります。

ここでフーリエ変換を思い出します。Rやmatlabで使われる式とは周期と積分範囲が異なるので注意してください。

フーリエ変換をかけます。公式により虚数とでの偏微分を無くす事ができます。

この式はスプリット・ステップ法（the split step method）で計算する事ができます。

オイラー法でを用いて書くと以下のようになります。

の形にして、両辺をで積分し、対数を消しての形に整理してから、を計算していますね。

上式の右辺第二項のフーリエ変換周りは注意してください。

] ]

フーリエ変換、逆フーリエ変換、フーリエ変換の三段構えです。

3. ソースコード

matlabのコードをRに翻訳した感じですが、意外に短いコードで計算する事ができます。計算速度は・・・遅いですね。なお精度はイマイチなので、ルンゲ＝クッタ法などに書き直した方が良いかも知れません。

# 補助的に使う関数
cosh <- function(x){ (exp(x)+exp(-x))/2 }
sech <- function(x){ 1/cosh(x) }
phi <- function(x, t, k=1){
  2*k^2*sech(k*(x-4*k^2*t))^2
}
ifft <- function(x){ fft(x, inverse=TRUE)/length(x) }


# 初期値
N <- 256
x <- seq(-10, 10, length.out=N)


# 波の初期値
# ソリトン波を指定しているが、他のモノでも良い
u1 <- phi(x, -1)


# 変化幅
delta_x <- x[2] - x[1]
delta_t <- 10/N^2 # エラーが出る場合は0.4/N^2で
# フーリエ変換の定義による周期の違いを2*piで調整
delta_k <- 2*pi/(N*delta_x)


# 離散フーリエ変換ではkの値が不連続になるので注意
# 参考: http://www.made.gifu-u.ac.jp/~tanaka/LectureNote/numerical_analysis.pdf P.109
k <- c(seq(0, N/2*delta_k, delta_k), seq(-(N/2-1)*delta_k, -delta_k, delta_k))


# t=1までnmax回波を動かす
tmax <- 1
nmax <- round(tmax/delta_t)
U = fft(u1)
for(n in 1:nmax){
  # first we solve the linear part
  U = U*exp(1i*k^3*delta_t)
  # then we solve the non linear part
  U = U - 3i*k*delta_t*(fft(Re(ifft(U))^2))
}
# 変化後の波を得る
u2 <- Re(ifft(U))


# 元と変化後の波を描画する
xlim = c(0, N)
ylim = c(0, max(u1, u2))
plot(u1, xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="x", ylab="u(t, x)", main="KdV方程式のプロット")
par(new=T)
plot(u2, xlim=xlim, ylim=ylim, type="l", xlab="", ylab="")
legend("topleft", lty=c(3, 1), col=c("black", "black"), legend=c("t=0", "t=1"))

4. 計算結果

アニメーションさせる気がしないぐらい時間がかかるので、2時点だけで。

Rでグラフの軸ラベルを回転させる

計量分析やシミュレーションに関係の無いデータで、Rでグラフを描いている人は少数派だと思います。
データ入力を表計算ソフトなどで行い、そこからテキスト・ファイルで保存するのに手間がかかる*1と言うのもあると思いますが、慣れていないと上手くレイアウトするのが大変だと言うのもあると思います。
例えば、X軸ラベルを縦書きにすると、以下のような操作が必要で大変です*2。

1. 元のグラフ

まずは、ありがちな折れ線グラフを描いてみます。軸のラベルの数が少ないので見栄えが良くありません。

2. 90度回転させる

ラベルを90度回転させると、軸ラベルの数が増えますが、出所と軸ラベルが重なります。

変更前コード

# 表示位置を12ごとにする。df$DATEは文字列型
at <- c(1, seq(1, length(df$DATE), by=12), length(df$DATE))
# 表示位置の部分だけ、df$DATEから抽出
labels <- df$DATE[at]
# ラベルの最後は空欄にする
labels[length(labels)] = ""
# X軸を表示する
axis(1, labels=labels, at=at, pos=1)

変更後コード

at <- c(1, seq(1, length(df$DATE), by=12), length(df$DATE))
labels <- df$DATE[at]
labels[length(labels)] = ""
# ラベルなしで軸だけ表示
axis(1, labels=FALSE, at=at, pos=0)
# par("usr")はグラフの表示位置を示すベクターで、3番目は下側の軸の位置になります。
# srt=-90がマイナス90度回転で、pos=4が左よりを示します。
text(x=at, par("usr")[3], labels=labels, srt=-90, pos=4, xpd=TRUE)

軸表示関数だけではラベルの回転をさせる事ができません。

3. 下側の余白を空け、出所の位置を移動させる

出所と軸ラベルが重なりを解消しましょう。

変更前コード

# marオプションは底辺，左側，上側，右側の順に余白の大きさを行の高さ（mex）で指定
par(mar=c(5, 4, 4, 2))

mtext("出所) 日本銀行時系列統計データ検索サイトのデータより編集", side=1, line=2.5, adj=0)

変更後コード

par(mar=c(8, 4, 2, 1))

mtext("出所) 日本銀行時系列統計データ検索サイトのデータより編集", side=1, line=5.5, adj=0)

4. まとめ

データセットと全体のコードが無いと分かりづらい面もあるため、最終形のファイルを公開しておきます。
将来的には改良されるかも知れませんが、現状ではRでグラフを描くのは、それなり手間隙がかかる事が分かると思います。プログラマー的には、文字列の縦横サイズを取得するメソッドが無いので、きっちり描画位置を決めづらい所も不満に残ります（strwidthとstrheightがあると言うツッコミがあったので取り消し）。

*1:Excelファイルのまま読み込ませるにはパッケージが必要だし、CSVファイルなどでも数値データにカンマが混じっていたりするとas.numeric(gsub(",", "", 列名))とするはめになります。

*2:90度回転でよければ、axis(1, las=2)とすれば楽です。

動学マクロ経済学と言う名の非線形連立方程式を解いてみる

動学マクロ経済学に関する蓮見氏のレクチャーノートを拝読していて、もう少しRのソースコード部分に補足説明があったほうが分かりやすい気がしたので、ラムゼー・モデル（Ramsey–Cass–Koopmans model）に関して紹介してみたいと思います。

このモデルは連続関数を仮定して変分法で解くこともできるのですが、最近は離散問題に置き換えてコンピューターで数値解を求める方向*1で教えているようです。これは、RBC/DSGEへの拡張が容易である*2ためでしょう。

何はともあれ、こういうモデルを計算してみると、マクロ経済学を直感的に理解できるようになると思います。解き方としてはモデルの式を整理して非線形連立方程式の形に直した上で、始点と終点を整理し、ニュートン・ラフソン法などで解くだけなので、実際にやってみましょう。

1. モデル式 → 非線形連立方程式

まずはモデルを確認しておきましょう。詳しくは蓮見氏のレクチャーノートを見てください。

1.1. 家計

元々のモデルではあまり特定化された関数は置いていない気がしますが、家計の目的関数になる効用関数を以下のようにおきます。

期の消費がになるわけですが、割引率があるので、将来の事はあまり気にしません。
制約条件は、以下のようになります。

資本に利率を乗じたものが利子収入、が賃金収入です。つまり、消費と投資は収入以下になると言う制約です。
来期の資本の量は今期の資本と投資の合計から、減価償却を引いたものになります。

これを制約条件に代入しておきましょう。

さて、制約付最大化問題を解くために、ラグランジュアンを置きます。

家計が決定できるのは、消費量と投資量なので、それぞれで偏微分をして一階条件を出します。の偏微分も忘れずに。

1.2. 企業

企業の生産関数は以下のように、コブ＝ダグラス型として起きます。

人口が増えないとしてと置きます。完全競争を仮定して、利子率と賃金を求めます。利子率は資本の限界生産性、賃金は労働の限界生産物に等しくなるわけですね。

生産物は利払いと賃金に分割されるので、以下のように企業の生産物と家計の収入が等しいくなります。

1.3. 非線形連立方程式

家計の効用最大化の一階条件の三つの式を整理しましょう。と、を消して、方程式の形式にしています。

1.4. 端点の整理

実は最適経路を求める問題なので、始点と終点が無いと計算できません。始点は初期保有ストック量、終点は定常状態となります。
定常状態はt期とt+1期の消費と資本ストックが同一と言う事になります。つまり、、となる点です。定常点をと置くと、非線形連立方程式の二式から以下のように整理できます。

2. カリブレーション

ここからが本題で、蓮見氏のレクチャーノートのソースコードにコメントが無いので、勝手に補足していきます。
数値計算する式は上の二式ですが、パラメーターや初期値を定めないといけませんし、非線形連立方程式を解いてくれるパッケージの使い方も確認する必要があります。

2.1. パラメーター

と以外の変数はパラメーターになります。任意でいいのですが、、、、となっていました。

2.2. 初期値

マクロ経済学と言うか数値解析の問題ですが、とに初期値が必要です。が1から30まで計算するとすると、2×30=60変数を初期化する必要があります。レクチャーノートでは何故か定常点が入っているのですが、0より大きな値であれば何でも問題ないです。

2.3. 計算までのソースコード

ようやくコーディングできるだけの情報が揃ったので計算してみましょう。なおレクチャーノートでは消費と投資の行列を作って崩していたのですが、ベクトルを分割するようにしています。

#
# 非線形連立方程式を解くパッケージを読み込む
#
library(nleqslv)


#
# パラメーターの設定
#
alpha <- 0.3
beta <- 0.99
delta <- 0.25
At <- 1.0


#
# 定常状態の計算
# （Kは値で、Cは作図のために関数で戻す）
#
steadystate <- function(){
  sK <- ((1/beta+delta-1)/alpha/At)^(1/(alpha-1))
  list(K=sK, C=function(K=sK){
    At*K^alpha - delta*K
  })
}
SS <- steadystate()


#
# カリブレーションを行う関数
# maxT:計算する期間の長さ, K01:資本ストック初期値
#
calibrate <- function(maxT, K01){
  # レキシカルスコーピングでcalibrate内で定義された変数maxTとK01は、
  # parseXとobjfunからアクセスできる一方でグローバル変数にならない


  #
  # 計算結果から消費と資本ストックを分ける
  # X:計算結果
  #
  parseX <- function(X) {
    # Xは(C[1], C[2], ... , C[maxT], K[2], K[3], ... , K[maxT + 1])と言うベクトル
  
    # Xの前半部分が各期の消費
    # 最後に定常状態を付け加える
    C <- c(X[1:maxT], SS$C())
  
    # Xの後半部分が各期の投資
    # 最初に初期値を付け加える
    K <- c(K01, X[(maxT+1):(2*maxT)])
  
    # リスト構造で値を戻す
    list(C=C, K=K)
  }
  
  #
  # 右辺ゼロの非線形連立方程式を表す関数
  # X:各期の消費と資本ストックを表すベクトル
  #
  objfun <- function(X) {
    # 消費と資本ストックに分ける
    X <- parseX(X)


    # 戻り値のために空ベクトルを創る
    rC <- numeric(maxT)
    rK <- numeric(maxT)
  
    # 線形連立方程式に値を代入していく
    # 2*maxT本の連立方程式になる事に注意
    for(t in 1:maxT) {
      rC[t] <- X$C[t+1]/X$C[t] - beta*(alpha*At*X$K[t+1]^(alpha-1)-delta+1)
      rK[t] <- X$K[t+1] - At*X$K[t]^alpha - (1-delta)*X$K[t]+X$C[t]
    }
  
    # 戻り値を返す
    return(c(rC, rK))
  }


  # 初期値は基本的に任意の値で構わないが、連立方程式内に割り算があるため0は避ける
  Xinit <- c(rep(1, maxT), rep(1, maxT))


  # 非線形連立方程式を解く
  r <- nleqslv(Xinit, objfun)


  # 結果から消費と資本ストックを分ける
  parseX(r$x)
}

# 期間の長さ
maxT <- 30


# 初期資本ストックは定常状態の半分として計算
X1 <- calibrate(maxT, 0.5*SS$K)


# 初期資本ストックは定常状態の倍として計算
X2 <- calibrate(maxT, 2*SS$K)

講義などで発散する経路が必要な場合は、以下のように計算します。

#
# 非定常経路の計算
#
inertia <- function(maxT, K01, C02){
  # 戻り値のために空ベクトルを作る
  C <- numeric(maxT)
  K <- numeric(maxT)
  K[1] <- K01
  C[1] <- C02
  for(t in 2:maxT){
    # 資本の増加式
    K[t] <- At*K[t-1]^alpha + (1 - delta)*K[t-1] - C[t-1]
    # オイラー方程式に従って来期の消費を決定
    C[t] <- beta*(alpha*( K[t-1]^alpha + (1 - delta)*K[t-1] - C[t-1])^(alpha-1) - delta + 1)*C[t-1]
  }
  list(C=C, K=K)
}

# 消費過少経路を計算
NSUC <- inertia(maxT, 0.5*SS$K, X1$C[1] - 0.02)

# 消費過剰経路を計算
NSOC <- inertia(maxT, 0.5*SS$K, X1$C[1] + 0.02)

2.4. 計算結果を描画する

X1$CとX1$Kを眺めてもいいのですが、見栄えが悪いので描画してみましょう。

#
# 定常状態を描画
#
curve(SS$C(x), 0, 3, ylim = c(0.5, 1.1), xlab = expression(K[t]), ylab = expression(C[t]), main="ラムゼーモデルの位相図", lty="dotted")
lines(cbind(c(SS$K, SS$K),c(0, 100)), lty="dotted")
points(SS$K, SS$C(), pch=19)


text(2.5, 0.7, expression(paste(Delta,"K = 0",sep="")), pos=4)
text(1.20, 1.0, expression(paste(Delta,"C = 0",sep="")), pos=4)


#
# 遷移を矢印で描画
# 矢印が密集していると見づらいので、間引く
#
step <- 3
for(i in seq(1, 12, step)){
  arrows(X1$K[i], X1$C[i], X1$K[i+step], X1$C[i+step], length = 0.1)
  arrows(X2$K[i], X2$C[i], X2$K[i+step], X2$C[i+step], length = 0.1)
}
for(i in seq(13, maxT - 1, step)){
  segments(X1$K[i], X1$C[i], X1$K[i+step], X1$C[i+step])
  segments(X2$K[i], X2$C[i], X2$K[i+step], X2$C[i+step])
}

描画される図形は以下のような感じになります。

発散経路を計算した場合は、以下のように続けて描画してください。

#
# 発散経路を描画
# すぐに発散するので間引かない
#
step <- 1
for(i in seq(1, maxT, step)){
arrows(NSOC$K[i], NSOC$C[i], NSOC$K[i+step], NSOC$C[i+step], length = 0.1, lty=1, col="red")
arrows(NSUC$K[i], NSUC$C[i], NSUC$K[i+step], NSUC$C[i+step], length = 0.1, lty=1, col="blue")
}

3. ラムゼー・モデルから分かること

最近の動学マクロ経済学は、もっと複雑なモデルが組まれていますし、パラメーターも実データから推定する世界になっていますが、このラムゼー・モデルが一つの出発的になっているのは間違いないと思います。動学マクロ経済学が、家計が将来を見通して意思決定を行う事を前提としていることが分かります。

*1:カリブレーションと言います。

*2:学生の数学力が低くて済むと言う利点もあるとは思います。