R中的NLS和對數週期冪律（LPPL）

這是我迄今爲止在R中完成的最具挑戰性的事情，因爲nls和LPPL對我來說都是相當新的。R中的NLS和對數週期冪律（LPPL）

下面是我一直在使用的腳本的一部分。 df是一個由兩列組成的數據框，日期和Y是S & P 500的收盤價。我不確定它是否相關，但是日期從01-01-01開始到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從這篇文章中的代碼，Using nls in R to re-create research。

謝謝你的幫助。

每jlhoward的評論，df.rda可以在這裏下載：https://drive.google.com/file/d/0B4xAKSwsHiEBb2lvQWR6T3NzUjA/edit?usp=sharing

來源

2014-02-15 mks212

看起來很熟悉......你需要顯示你的數據幀'df'，或者你用來下載和處理它的代碼。否則，不可能提供有用的答案。警告可能是因爲xx> pars $ tc。 – jlhoward

已發佈。以前沒有信任你的道歉，已被糾正。 – mks212

首先，一對夫婦的小事情：

兩個nls(...)和nls.lm(...)需要的數值參數，沒有日期。所以你必須以某種方式轉換。我剛剛添加了days列，這是自數據開始以來的天數。
您對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-18 06:54:08 jlhoward

R中的NLS和對數週期冪律（LPPL）

回答

相關問題