餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

Rのグラフで余白の調整

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で多体問題を数値解析

R言語

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

1. 計算する式

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

m_i \frac{d^2r_i}{dt^2}=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}}

Gが万有引力定数、mが質量、rが距離、nは質点の数、iは質点、tが時点を表します。
力の成分をxとyに分解します。z成分は省略。

m_i \frac{d^2x_i}{dt^2}=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}} \cdot \frac{x_j-x_i}{r_{ij}}

m_i \frac{d^2y_i}{dt^2}=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}} \cdot \frac{y_j-y_i}{r_{ij}}

r_{ij} = \sqrt{(x_j-x_i)^2 + (y_j-y_i)^2}

両辺をm_iで割って、x2=\frac{d^2x_i}{dt^2}、y2=\frac{d^2y_i}{dt^2}と置いて、常微分方程式に変形しましょう。

x_1=\frac{dx_i}{dt}
y_1=\frac{dy_i}{dt}
x_2=-\sum^n_{j=1,j \neq i}G\frac{m_j}{r^2_{ij}} \cdot \frac{x_j-x_i}{r_{ij}}
y_2=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}} \cdot \frac{y_j-y_i}{r_{ij}}

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

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方程式を数値解析でプロットしてみる

R言語

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

1. KdV方程式の形状

まずは式の確認です。

\frac{\partial}{\partial t} u(t,x) + 6u(t,x)\frac{\partial}{\partial x}u(t,x) + \frac{\partial^3}{\partial x^3} u(t,x) =0

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

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

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

\frac{\partial}{\partial t} u(t,x) + 3\frac{\partial}{\partial x} u(t,x)^2 + \frac{\partial^3}{\partial x^3} u(t,x) = 0

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

F\bigg(u(k) \bigg ) = \hat{u(k)} = \int^{\infty}_{-\infty} u(x) e^{-ikx} dx

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

\frac{\partial}{\partial t} \hat{ u(k) } + 3ik \hat{ u(k)^2 } - ik^3 \hat{ u(k) } = 0

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

\frac{\partial}{\partial t} \hat{ u(k) } = - 3ik \hat{ u(k)^2 } + ik^3 \hat{ u(k) }

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

\hat{ u_1(k, t + \Delta t) } = \hat{ u(k, t) } e^{ik^3 \Delta t}

\partial \hat{u}/\partial t \cdot 1/\hat{u}=ik^3の形にして、両辺をtで積分し、対数を消してu(k, t)=C e^{ik^3 t}の形に整理してから、u(k, t + \delta t)を計算していますね。

\hat{ u(k, t + \Delta t) } = \hat{ u_1(k, t + \Delta t) } - 3ik \Delta t \hat{ u_1(k, t + \Delta t)^2 }

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

\hat{ u(k, t + \Delta t) } = \hat{ u_1(k, t + \Delta t) } - 3ik \Delta t F [ F^{-1} [ \hat{ u_1(k, t + \Delta t)^2 } ] ]

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

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言語

計量分析やシミュレーションに関係の無いデータで、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言語

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

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

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

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

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

1.1. 家計

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

\max_{C} U_t = \sum^{\infty}_{i=t} \beta^{i-t} \ln(C_{i})

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

\mbox{s.t. } r_t K_i + w_i - C_i - I_{i} \geq 0

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

K_{i+1}=K_i - \delta K_i + I_i

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

\mbox{s.t. } r_tK_i + w_i + K_i - \delta K_i - C_i - K_{i+1} \geq 0

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

\mathit{L}=\sum^{\infty}_{i=t} \beta^{i-t} \bigg[ \ln(C_{i}) + \lambda_i\{ r_tK_i + w_i + K_i - \delta K_i - C_i - K_{i+1} \} \bigg]

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

\frac{1}{C_i}-\lambda_i=0

\beta(r_{t+1}+1-\delta)\lambda_{i+1} - \lambda_i=0

r_tK_i + w_i + K_i - \delta K_i - C_i - K_{i+1}=0

1.2. 企業

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

Y_t=A_t K^\alpha L^{1-\alpha}

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

r_t=\alpha A_t K_t^{\alpha-1}

w_t= A_t K_t^\alpha-r_tK_t = (1-\alpha) A_t K_t^{\alpha}

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

Y_t=A_t K_t^\alpha=r_tK_i + w_i

1.3. 非線形連立方程式

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

\frac{C_{t+1}}{C_t}-\beta(\alpha A_{t+1} K_{t+1}^{\alpha-1}-\delta+1)=0

K_{t+1}-(A_t K_t^\alpha + (1-\delta)K_t - C_t)=0

1.4. 端点の整理

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

K^*=\bigg [ \frac{\frac{1}{\beta} + \delta - 1}{\alpha A_t} \bigg ]^{\frac{1}{\alpha - 1}}
C^*= A_t \cdot K^{*\alpha} - \delta \cdot K^*

2. カリブレーション

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

2.1. パラメーター

C_tとK_t以外の変数はパラメーターになります。任意でいいのですが、\alpha=0.3、\beta=0.99、\delta=0.25、A_t=1となっていました。

2.2. 初期値

マクロ経済学と言うか数値解析の問題ですが、C_tとK_tに初期値が必要です。tが1から30まで計算するとすると、2×30=60変数を初期化する必要があります。レクチャーノートでは何故か定常点(C^*, K^*)が入っているのですが、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:学生の数学力が低くて済むと言う利点もあるとは思います。

Rのグラフにスプライン曲線を描きこむ

R言語

TTFなどのフォントもそうですし、CADなどグラフィックス分野の方が使う事が多いと思いますが、間欠的なデータから連続した線を描画したいときは実用上は多々あります。ラグランジュ補間、ニュートン補間などバリエーションは幾つもありますが、Rではスプライン補間をする事が容易です。
「R でグラフ中にプロットされてない値を予測してみる」を見ていたら簡単そうなので、試してみました。y=\frac{\sin(x)}{x}と言う簡単そうで積分するとトリッキーで大変なルベーグ積分不可の関数の極値をつないでみます。
まずは、y=\frac{\sin(x)}{x}を描いてみます。

次に、極値を求めます。

最後に、極大値と極小値をスプライン曲線で結びます。

意味が良く分からないものの、それらしいグラフになりました。ソースコードは以下になります。

# 関数を定義
f <- function(x){
  sin(x)/x
}


# fの一階微分
df <- function(x){
  (cos(x)*x-sin(x))/x^2
}


# fの二階微分
ddf <- function(x){
  (-sin(x)*x*2*x - (cos(x)*x-sin(x))*2*x)/x^4
}


# xの変分
dx <- 0.01
# xを生成（始点を僅かに正にすることでyが∞を回避）
x <- seq(1e-6, 30, dx)
# yを計算
y <- f(x)


getExtremum <- function(x){
  dy <- df(x)
  n <- length(dy)
  # 微分係数dy[i]*dy[i+1]がゼロ以下の場合は符号が逆転しており極値近傍と見なせる
  x[dy[-n]*dy[-1] <= 0]
}
extremum <- getExtremum(x)


# ddf(x[i])が負ならば、極大値
ulim <- c(extremum[ddf(extremum)<0])
# ddf(x[i])が正ならば、極小値
llim <- c(extremum[ddf(extremum)>0])


# スプライン関数を導出
sp_u <- smooth.spline(ulim, f(ulim))
sp_l <- smooth.spline(llim, f(llim))


# xとyのそれぞれ上下限
ylim <- c(-0.5, 1)
xlim <- c(min(x), max(x))


# 余白位置などを調整
par(mgp=c(1, 1, 0), mai=c(1,1,1,1), mar=c(4,3,3,1))


# 元データをプロット
curve(f(x), xlim=xlim, ylim=ylim, type="l", xlab="x", ylab=expression(sin(x)/x), axes=FALSE, main="スプライン関数のお試し")


# スプラインを描画／predict()関数で補間点を計算
par(new=T); plot(predict(sp_u, x), xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="", ylab="", axes=FALSE, col="red2")
par(new=T); plot(predict(sp_l, x), xlim=xlim, ylim=ylim, type="l", lty="dotted", xlab="", ylab="", axes=FALSE, col="red2")


# スプライン関数の元もプロットしておく
par(new=T); plot(ulim, f(ulim), xlim=xlim, ylim=ylim, type="p", lty="dotted", xlab="", ylab="", axes=FALSE, bg="gray")
par(new=T); plot(llim, f(llim), xlim=xlim, ylim=ylim, type="p", lty="dotted", xlab="", ylab="", axes=FALSE, bg="gray")


# 軸を描画
axis(1, at=seq(0, 30, 10), pos=ylim[1])
axis(2, at=c(seq(ylim[1], ylim[2], 0.25)), pos=0)

極大値と極小値を計算するほうが手間がかかっていますね。

Rで偏微分方程式

R言語

速度的に実用性が低い感じもするのですが、Rで偏微分方程式の数値解を計算することもできます。微分方程式一般に言えることですが、せっせと計算して閉じた解を求めるのは苦労するので、数値解は実用的には重要です。
deSolveパッケージを使えば三元連立式ぐらいまでは容易で、使い方はサンプル・コードを見ればよいのですが、恐らく最も見慣れている偏微分方程式である波動方程式のサンプルを書いてみました。
deSolveの旧バージョンであるReacTranパッケージのリファレンスに例としては上がっているわけですが、パッケージの構成が変わってそのままでは動かないので、もしかしたら参考になるかも知れません。

1. 計算の概要

まずは波動方程式を眺めてみましょう。波動関数はu(t,x)で、xを位置、tを時間とします。cは波動の速度。
\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}
偏微分方程式のままだと解けないので、常微分方程式に変形します。
\frac{\partial u_1}{\partial t} = u_2
\frac{\partial u_2}{\partial t} = c^2 \frac{\partial^2 u_1}{\partial x^2}
初期パラメーターを以下のようにします。力学的振動のときと異なって、ベクトルで大量の値がセットされることに注意してください。
u_1(0,x)=exp{-0.05x^2}
u_2(0,x)=u(t,\infty)=u(t,-\infty)=0
c=1
この式からそのままコーディングできる人は、かなり数値演算に慣れていると思います。

2. 実際の演算コード

以下が実際のコードになります。ポイントは恐らく、二階微分を行うところになると思います。数値微分のちょっとした知識が要ります。deSolveパッケージの癖で、複数の変数のベクトルが、一つのベクトルにまとめられているのも、注意しないと混乱するかも知れません。

# deSolveパッケージを読み出す
library("deSolve")
# 位置xの刻み幅は0.2
dx <- 0.2
# xの最小値と最大値
xlim <- c(-70, 70)
# -40から40まで位置xのベクトルを得る
x <- seq(xlim[1], xlim[2], dx)
# xの数を保存しておく
N <- length(x)
# 波の高さの初期値
uini <- exp(-0.05 * x^2)
# 波の高さの変化分の初期値
vini <- rep(0, N)
# パラメーターは同じベクターにしておく
yini <- c(uini, vini)
# 時間tのベクトルを得る
times <- seq (from = 0, to = 50, by = 1)
# 常微分方程式の形式で、波動方程式を表現する
wave <- function(t, y, parms, dx, N) {
 # ベクターを二つに分ける
 u1 <- y[1:N]
 u2 <- y[(N+1):(2*N)]
 # 波の高さの変化を計算
 du1 <- u2
 # 波の高さの変化分の変化の計算
 # 二階微分を行う。両端の変化率は境界条件によりゼロ。
 du2 <- c(0, (diff(u1[-1],1) - diff(u1[-N],1))/dx^2, 0)
 # 変化を戻り値で返す
 return(list(c(du1, du2)))
}
# ode.1D関数は一次元の微分方程式を扱う。
# nspecは計算する変数の数で、ここではu1とu2で二つ。
# 計算手法はデフォルトだと収束しなくなるので、ode45
out <- ode.1D(yini, times, wave, parms = 0, nspec = 2, dx=dx, N=N, method="ode45", atol=1e-14, rtol=1e-14)

ちょっと時間がかかると思いますが、気長に待てば計算できるでしょう。

3. 演算結果をプロットしてみる

出てきた計算結果は以下のようにプロットしたり、中身を見ることができます。

# 波の変化をプロットしてみる
matplot.1D(out, which=1, subset = time %in% seq(0, 50, by = 10), type = "l", col = c("black", rep("darkgrey", 5)), lwd = 2, grid = x, xlim = xlim, main="u(x, t)")
legend("topright", lty = 1:6, lwd = 2, col = c("black", rep("darkgrey", 5)), legend = paste("t = ",seq(0, 50, by = 10)))

# 等高線プロットをしてみる
image(out, xlab = "t", ylab = "x", main = "PDE", add.contour = FALSE)

左がu(t, x)、右がduになります。色が波動の高さで、赤色が高く青色が低くなります。表題などを綺麗にする方法は調査中。

# t=30のときの大きさを波の形状を見てみる。ylim重要。
plot(out[30,][1:N], ylim=c(0, 1), xlab="x", ylab="u")
# 中央地点x=2*N/3の波の高さの変化を見る
plot(out[,2*N/3][1:length(times)], ylim=c(0, 1), xlab="t", ylab="u")

コンピューターを使っている気になれますね。