0

見積もりの​​問題で立ち往生しているのではないかと心配しています。

X と Y の 2 つの変数があります。Y は、X の n 個のラグ値の加重和によって説明されます。私の目的は、次の 2 つのパラメーター c(alpha0,alpha1) を推定することです。

Yt = ( ( alpha0 + alpha1 * j ) * Xt-j ) の j=1 から n までの合計

ここに画像の説明を入力

ここで、Xt-j は X の j 番目の遅れを表します。

このアプローチを思いついたのは、追加された X のラグごとに 1 つのパラメーターを推定するのではなく、重みの勾配を推定する方がよいと考えたからです (n を非常に大きく設定するつもりです)。

モデルにノイズ ut が追加されます。これは、平均ゼロと標準偏差シグマで正規分布していると想定されます。

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 = ( ( alpha0 + alpha1 * j ) * Xt-j ) の j=1 から n までの合計

= (アルファ+ベータ*1)*Xt-1 + (アルファ+ベータ*2)*Xt-2 + (アルファ+ベータ*3)*Xt-3 + ... + (アルファ+ベータ*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関数は、ラグ関数が行うべきだと思っていたことに役立ちます。

> 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 には不完全なデータがあるため、期間 1 ～ 6 で何をすべきかがあいまいになっているという事実を受け入れて計画する必要があります。磨耗現象に特定の形状を適用しない限り、2つ以上のラグ期間がある場合に2つのパラメーターのみを推定する方法をあなたの式から理解することはできません...おそらく線形...あなたは言いませんでした。

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 に答える