1

这是迄今为止我在 R 中所做的最具挑战性的事情，因为 nls 和 LPPL 对我来说都是相当新的。

下面是我一直在使用的脚本的一部分。df 是一个数据框，由两列组成，日期和 Y，它们是标准普尔 500 指数的收盘价。我不确定它是否相关，但日期从 01-01-2003 到 12-31-2007。

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 屏幕上找到，http://www.chronostraders.com/wp-content/uploads/2013/08/Research_on_LPPL.pdf

你知道我哪里错了吗？这适用于不同的模型，我更改了新方程的代码。这篇文章中的这段代码归功于 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 方程与 Eqn 不同。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和omeega具有p < 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 回答