2013-04-01 2 views
0

Я боюсь, что я застрял в задаче оценки.Оценка параметров в R из взвешенной суммы отстающих переменных

У меня есть две переменные, X и Y. У объясняется взвешенной суммы п отставал значения X. Моя цель состоит в том, чтобы оценить два параметра с (alpha0, альфа1) в:

Yt = СУММ от J = 1 до п ((alpha0 + альфа1 * J) * Х-J)

enter image description here

где Xt-J обозначает запаздывание его из X.

я пришел с этим подходом потому что я думал, что было бы неплохо оценить наклон весов вместо того, чтобы оценивать один параметр для каждого запаздывания X (я намерен установить n очень большой).

К модели шума ut добавлен, который считается нормально распределенным со средним значением нуля и стандартным отклонением сигма.

Предполагая, что я хотел бы установить n = 510, тогда мне нужна оригинальная серия и 510 отстающих серий. Чтобы избежать любых НС в серии, я преобразовываю исходные данные в «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 = Сумма от у = 1 до п ((alpha0 + альфа1 * J) * Х-J)

= (альфа + бета * 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 

Может кто-нибудь, пожалуйста, помогите мне в этом?

+0

Похоже, вам нужно подогнать некоторую модель «ARIMA» или, скорее, модель «MA» ... – agstudy

+0

Возможно, вам понадобится некоторый шум в вашей модели, 'y [t] = ... + epsilon [t]'. Ваша логарифмическая вероятность должна использовать данные 'y' и, если шум гауссов, он должен содержать сумму квадратов. Избегайте функции 'lag': (если вы не являетесь« временными рядами »xts или' zoo'), он никогда не делает то, что вы хотите. В вашей функции 'logl'' 'NULL': вы вычисляете что-то внутри цикла, , но вы ничего не делаете с результатом, поэтому он отбрасывается. –

+0

Привет Винсент! Спасибо вам за ваш ответ. Я действительно занимаюсь данными ts здесь, поэтому, я думаю, я могу держать лагами? Поэтому вы попросили меня включить в функцию остатки. Является ли ССБ фактически тем, что должно быть минимизировано с помощью оптимизации? – Chris437

ответ

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, поскольку они будут иметь неполные данные. Я не могу понять из вашей формулы, как можно оценить только два параметра, если у вас более двух периодов задержки, если вы не примените определенную форму к феномену износа .... линейный, возможно ... вы не сказали.

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 и факторы, вплоть до 1/7-й силы для X_ (трет-7). Это немного отличается от того, что выражала ваша формула, поскольку у нее не было x_t covariate, но вы должны иметь возможность построить «наклон» из оцененных коэффициентов.