1

これは、nls と LPPL の両方が私にとってかなり新しいものであるという点で、これまで R で行った中で最も困難なことです。

以下は、私が取り組んできたスクリプトの一部です。df は、S&P 500 の終値である日付と Y の 2 つの列で構成されるデータ フレームです。関連性があるかどうかはわかりませんが、日付は 2003 年 1 月 1 日から 2007 年 12 月 31 日までです。

f <- function(pars, xx) {pars$a + pars$b*(pars$tc - xx)^pars$m * 
                    (1 + pars$c * cos(pars$omega*log(pars$tc - xx) + pars$phi))} 
# residual function
resids <- function(p, observed, xx) {df$Y - f(p,xx)}
# fit using Levenberg-Marquardt algorithm
nls.out <- nls.lm(par=list(a=1,b=-1,tc=5000, m=0.5, omega=1, phi=1, c=1 ), fn = resids, 
              observed = df$Y, xx = df$days)
# use output of L-M algorithm as starting estimates in nls(...)
par <- nls.out$par

nls.final <- nls(Y~a+b*(tc-days)^m * (1 + c * cos(omega * log(tc-days) + phi)),data=df, 
             start=c(a=par$a, b=par$b, tc=par$tc, m=par$m, omega=par$omega, phi=par$phi,         c=par$c))
summary(nls.final) # display statistics of the fit 
# append fitted values to df
df$pred <- predict(nls.final)

実行すると、次のメッセージが表示されます。

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates
In addition: Warning messages:
1: In log(pars$tc - xx) : NaNs produced
2: In log(pars$tc - xx) : NaNs produced

LPPL の公式は、この PDF ファイルの 5 番目の画面にあります。

私がどこで間違っているか知っていますか？これは別のモデルでは正しく機能していたので、新しい方程式のコードを変更しました。この投稿のこのコードについては、jlhoward の功績によるもので、R で nls を使用して研究を再作成します。

ご協力ありがとうございました。

jlhoward のコメントによると、df.rda はここからダウンロードできます: https://drive.google.com/file/d/0B4xAKSwsHiEBb2lvQWR6T3NzUjA/edit?usp=sharing

4

1 に答える 1

3

まず、いくつかのマイナーなことを説明します。

  1. との両方nls(...)に、nls.lm(...)日付ではなく数値引数が必要です。したがって、何らかの方法で変換する必要があります。daysデータの開始からの日数である列を追加しました。
  2. F の式は、式とは異なります。参考に1があったので、合わせるように変更しました。

*

f <- function(pars, xx) 
         with(pars,(a + (tc - xx)^m * (b + c * cos(omega*log(tc - xx) + phi))))

ここで大きな問題があります: あなたの最初の見積もりは、LM 回帰が収束しないようなものです。その結果、 の値nls.out$parは安定した推定値ではありません。これらを の開始点として使用すると、nls(...)同様に失敗します。

nls.out <- nls.lm(par=list(a=1,b=-1,tc=5000, m=0.5, omega=1, phi=1, c=1 ),
                  fn = resids, observed = df$Y, xx = df$days)
# Warning messages:
# 1: In log(pars$tc - xx) : NaNs produced
# 2: In log(pars$tc - xx) : NaNs produced
# ...
# 7: In nls.lm(par = list(a = 1, b = -1, tc = 5000, m = 0.5, omega = 1,  :
#   lmdif: info = -1. Number of iterations has reached `maxiter' == 50.

一般に、何が起こったのかを調べて確認する必要がnls.out$statusあります。nls.out$message

7 つのパラメーターを持つ複雑なモデルがあります。残念ながら、これは回帰が多くの極小値を持つ状況につながります。したがって、収束につながる見積もりを提供しても、「役に立たない」可能性があります。検討：

nls.out <- nls.lm(par=list(a=1,b=1,tc=2000, m=-1, omega=1, phi=1, c=1 ), 
                  fn = resids, observed = df$Y, xx = df$days, 
                  control=nls.lm.control(maxiter=10000, ftol=1e-6, maxfev=1e6))
par <- nls.out$par
par
plot(df$Date,df$Y,type="l")
lines(df$Date,f(par,df$days))

これは安定した結果 (極小値) ですが、c比較すると非常に小さいためb、振動は目に見えません。一方、これらの開始推定値は、基準とかなり厳密に一致する適合を生成します。

nls.out <- nls.lm(par=list(a=0,b=1000,tc=2000, m=-1, omega=10, phi=1, c=200 ), 
                  fn = resids, observed = df$Y, xx = df$days, 
                  control=nls.lm.control(maxiter=10000, ftol=1e-6, maxfev=1e6))

これにより、 との収束につながるパラメーター推定値が生成されますnls(...)が、要約では、パラメーターの推定が不十分であることが示されています (とtcのみ)。omeegap < 0.05

nls.final <- nls(Y~a+(tc-days)^m * (b + c * cos(omega * log(tc-days) + phi)),
                 data=df, start=par, algorithm="plinear",
                 control=nls.control(maxiter=1000, minFactor=1e-8))
summary(nls.final)

最後に、参考文献 (大恐慌ではなく大恐慌をモデル化していることは確かです) に非常に近い開始推定値を使用すると、さらに優れた結果が得られます。

nls.out <- nls.lm(par=list(a=600,b=-266,tc=3000, m=.5,omega=7.8,phi=-4,c=-14), 
                  fn = resids, observed = df$Y, xx = df$days, 
                  control=nls.lm.control(maxiter=10000, ftol=1e-6, maxfev=1e6))

于 2014-02-18T06:54:08.173 に答える