这是迄今为止我在 R 中做过的最具挑战性的事情，因为 nls 和 LPPL 对我来说都是相当新的。
以下是我一直在使用的脚本的一部分。 df 是一个由两列组成的数据框，Date 和 Y，它们是 S&P 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 文件的第五屏中找到，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
首先，有一些小事情：
nls(...) and nls.lm(...)需要数字参数，而不是日期。所以你必须以某种方式进行转换。我刚刚添加了一个days列是自数据开始以来的天数。*
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 and 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 and omeega have 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))
