Vasicek Modelで短期金利パスを発生させる

http://www.amazon.com/dp/0387004513/のP１１０をそのまま実装。

#create short rate path by one factor stationary Vasicek model
#alpha:緩和時間の逆数
#b:回帰水準
#dt:時間刻み幅
#endTime:終了時間
#initialRate:初期時点での金利
StationaryVasicek <- function(alpha,b,sigma,dt,endTime,initialRate)
{
    shortRate <- c(initialRate,numeric(endTime/dt))
    randNormal <- rnorm(endTime/dt)
    for(i in 1:(endTime/dt)){
        shortRate[i+1] <- exp(-alpha*dt)*shortRate[i] + b*(1-exp(-alpha*dt)) + 
                            sigma*sqrt(1/(2*alpha)*(1-exp(-2*alpha*dt)))*randNormal[i]
    }
    return(shortRate)
}

金利パスを発生させて描画してみる

shortRate <- StationaryVasicek(alpha=0.01,b=1,sigma=1,dt=0.01,endTime=10,initialRate=1)
plot(shortRate)

と、

な感じで短期金利の変動がグラフ化される。
forループを回さない設計がしたいです。

〜追記（2012/5/20）〜
さすがに初稿をUPしてから2年近くも経てば私も少しは進歩しているようで一応forループなし版が出来ました。こちゃこちゃしているのであまりオススメできる書き方ではないですが、ご参考に。

#１ステップ更新するための関数定義
NextStep <- function(alpha, b , sigma, dt){
  function(r, rand){
    exp(-alpha*dt)*r + b*(1-exp(-alpha*dt)) + 
      sigma*sqrt(1/(2*alpha)*(1-exp(-2*alpha*dt)))*rand
  }
}
#金利モデルのパラメーターを束縛＆１ステップ計算用の関数オブジェクト生成
next.step <- NextStep(alpha = 0.01, b = 1, sigma = 1, dt = 1.0)
#パスの生成
Reduce(function(x, y){next.step(x, y)}, x = rnorm(10), init = 1, accumulate = TRUE)

最後のパス生成の結果は当然上の関数（StationaryVasicek関数）と一致していて、乱数のシードを固定する事で以下のようにチェックできる。

> set.seed(10)
> StationaryVasicek(alpha=0.01,b=1,sigma=1,dt=1.0,endTime=10,initialRate=1)
 [1]  1.0000000  1.0186528  0.8351321 -0.5277298 -1.1087130 -0.7946524 -0.3889419 -1.5771826 -1.9134044
[10] -3.5029886 -3.7133844
> 
> set.seed(10)
> Reduce(function(x, y){next.step(x, y)}, x = rnorm(10), init = 1, accumulate = TRUE)
 [1]  1.0000000  1.0186528  0.8351321 -0.5277298 -1.1087130 -0.7946524 -0.3889419 -1.5771826 -1.9134044
[10] -3.5029886 -3.7133844

バッチリ一致している。