0

恐怕我遇到了估计问题。

我有两个变量，X 和 Y。Y 由 X 的 n 个滞后值的加权和来解释。我的目标是估计以下两个参数 c(alpha0,alpha1)：

Yt = 从 j=1 到 n 的总和 ( ( alpha0 + alpha1 * j ) * Xt-j )

在此处输入图像描述

其中 Xt-j 表示 X 的第 j 个滞后。

我提出了这种方法，因为我认为估计权重的斜率而不是为每个添加的 X 滞后估计一个参数是一个好主意（我打算将 n 设置得非常大）。

模型噪声 ut 被添加，假设它是正态分布的，均值为 0，标准差为 sigma。

假设我想设置n=510，那么我需要原始系列和 510 滞后系列。为了避免系列中的任何 NA，我将原始数据转换为“data_chopped”，仅包含前 510 个观察值被删除后的观察值，以及矩阵“data_lagged”，其中每列代表一个滞后序列：

library(stats)
data<-arima.sim(n=10000,list(ar=0.15,ma=0.1),mean=0.5)

data_chopped<-data[511:length(data)]

data_lagged<-matrix(nrow=length(data_chopped),ncol=510)
for (i in 1:510){
data_lagged[,i]<-head(data,-i)[(511-i):length(head(data,-i))]
}

#Check result:
cbind(data_chopped,data_lagged[,1:3])
#data_lagged[,1] is the first lag of the original data, data_lagged[,2] is the second lag, and so on. No NAs whatsoever to deal with

为了演示我的对数似然函数和生成的系列的“工作顺序”，我首先想拟合一个 AR(3) 模型：

logl<-function(sigma,alpha,beta,gamma){
-sum(log((1/(sqrt(2*pi)*sigma)) * exp(-((
data_chopped
-alpha*data_lagged[,1]
-beta*data_lagged[,2]
-gamma*data_lagged[,3]
)^2)/(2*sigma^2))))
}

library(stats4)
mle(logl,start=list(sigma=1,alpha=0,beta=0,gamma=0),method="L-BFGS-B")

当我现在尝试以同样的方式估计我的模型时，它就不起作用了。我从来没有在对数似然函数中使用循环，这就是我刚刚写出上述模型的原因。所以，

Yt = 从 j=1 到 n 的总和 ( ( alpha0 + alpha1 * j ) * Xt-j )

= (alpha+beta*1)*Xt-1 + (alpha+beta*2)*Xt-2 + (alpha+beta*3)*Xt-3 + ... + (alpha+beta*510)*Xt -510

logl<-function(sigma,alpha,beta){
-sum(log((1/(sqrt(2*pi)*sigma)) * exp(-((
data_chopped
-(alpha + beta*1)*data_lagged[,1]
-(alpha + beta*2)*data_lagged[,2]
-(alpha + beta*3)*data_lagged[,3]
-(alpha + beta*4)*data_lagged[,4]
-(alpha + beta*5)*data_lagged[,5]
...
-(alpha + beta*510)*data_lagged[,510]
)^2)/(2*sigma^2))))
}

library(stats4)
mle(logl,start=list(sigma=1,alpha=0.5,beta=0),method="L-BFGS-B")
Error in optim(start, f, method = method, hessian = TRUE, ...) : 
L-BFGS-B needs finite values of 'fn'

如果我只尝试几行，我不会收到错误：

logl<-function(sigma,alpha,beta){
-sum(log((1/(sqrt(2*pi)*sigma)) * exp(-((
data_chopped
-(alpha + beta*1)*data_lagged[,1]
-(alpha + beta*2)*data_lagged[,2]
-(alpha + beta*3)*data_lagged[,3]
-(alpha + beta*4)*data_lagged[,4]
-(alpha + beta*5)*data_lagged[,5]
)^2)/(2*sigma^2))))
}

library(stats4)
mle(logl,start=list(sigma=1,alpha=0.5,beta=0),method="L-BFGS-B")
Call:
mle(minuslogl = logl, start = list(sigma = 1, alpha = 0.5, beta = 0), 
method = "L-BFGS-B")

Coefficients:
 sigma      alpha       beta 
1.07797708  0.26178848 -0.04378526

有人可以帮我解决这个问题吗？

4

1 回答 1

1

我会回应不使用该lag功能的建议。它的作者和早期用户可能知道它的作用，但我们其他人在它没有达到预期的情况下有过糟糕的体验。我发现这个embed函数对我认为 lag 函数应该做的事情很有用。

> embed(1:8, 3)
     [,1] [,2] [,3]
[1,]    3    2    1
[2,]    4    3    2
[3,]    5    4    3
[4,]    6    5    4
[5,]    7    6    5
[6,]    8    7    6

假设您想回顾当前时间之前的 6 次并按行进行计算。您需要接受并计划这样一个事实，即现在应该对第 1-6 期做什么是模棱两可的，因为它们将有不完整的数据。我无法从您的公式中弄清楚，当您有两个以上的滞后期时，如何仅估计两个参数，除非您将某些特定形状应用于磨损现象......也许是线性的......你没有说。

dfrm <- data.frame(y=rnorm(20), x=rnorm(20) )
dfrm$embx<- matrix(NA, ncol=7, nrow=20)
dfrm$embx[7:20, ] <- embed(dfrm$x, 7) * rep( (7:1)/7, each=14)
lm( y[7:20] ~ embx[7:20,], data=dfrm )

Call:
lm(formula = y[7:20] ~ embx[7:20, ], data = dfrm)

Coefficients:
  (Intercept)  embx[7:20, ]1  embx[7:20, ]2  embx[7:20, ]3  embx[7:20, ]4  embx[7:20, ]5  
       0.3065        -0.2371         0.9504         0.8601         0.5484         0.6621  
embx[7:20, ]6  embx[7:20, ]7  
       1.1619         4.8338

这使用“全强度” x_t 并将 x_(t-7) 的强度降低到 1/7。这与你的公式表达的有点不同，因为它没有 x_t 协变量，但你应该能够从估计的系数中构造一个“斜率”。

于 2013-04-02T00:43:39.783 回答