Модель многомерного пространства состояний закон Олунса

Question

Я пытаюсь оценить уравнение закона Окуня с помощью dlm, используя пакет dlm в R. Я могу оценить не изменяющуюся во времени модель, используя nls, следующим образом:

const_coef <- nls(formula = dur~ b1*dur_lag1 + b2*(d2lgdp-b0) + b3*d2lrulc_lag2 ,
start = list(b0 =0.1, b1=0.1, b2=0.1, b3=0.1),
data = mod_data) 

модель DLM я хочу быть в состоянии оценить позволяет b1 а также b0 в вышесказанном следовать случайным прогулкам. Я могу сделать это в Eviews, объявив уравнение измерения и добавив состояния (ниже приведен код, предоставленный авторами оригинальной статьи, который я могу воспроизвести:

'==========================
' SPECIFY THE KALMAN FILTER
'==========================

'Priors on state variables
vector(2) mprior
    mprior(1) = 4           'Prior on starting value for trend GDP growth (annual average GDP growth over 1950s)
    mprior(2) = 0           'Prior on starting value for lagged dependent variable
sym(2) vprior
    vprior(1,1) = 5         'Prior on variance of trend GDP growth   (variance of annual GDP growth over 1950s)
    vprior(2,2) = 1         'Prior on variance of lagged dependent variable

'Specify coefficient vector
coef(8) ckf

'Declare state space
sspace ss1
ss1.append dur = lag*dur(-1) + ckf(2)*(d2lgdp-trend)+ckf(3)*D2LRULC(-2)+[var=exp(ckf(4))]   'Measurement equation
ss1.append @state trend = 1*trend(-1) + [var = exp(ckf(5))]                                                 'State equation for trend GDP growth (random walk)
ss1.append @state lag = 1*lag(-1) + [var = exp(ckf(6))]                                                     'State equation for lagged dependent variable (random walk)

'Apply priors to state space
ss1.append @mprior mprior
ss1.append @vprior vprior

'Set parameter starting values
param ckf(2) -0.0495 ckf(3) 0.01942 ckf(4) -2.8913 ckf(5) -4.1757 ckf(6) -6.2466        'starting values for parameters

'=====================
' ESTIMATE THE MODEL
'=====================

'Estimate state space
smpl %estsd %ested          'Estimation sample
ss1.ml(m=500,showopts)      'Estimate Kalman filter by maximum likelihood
freeze(mytab) ss1.stats

Я действительно не уверен, как это сделать с пакетом dlm. Я пробовал следующее:

buildSS <- function(v){


  dV <- exp(v[1])               # Variance of the measurment equation (ckf4)
  dW <- c(exp(v[2]),            # variance of the lagged dep  (ckf6)
          0,                    # variance of the coef on d2lgdp ckf(2) set to 0
          0,                    # variance of the coef on d2lrulc ckf(3) set to 0
          exp(v[3])             # variance of the random walk intercept (ckf5)
           )

  beta.vec <- c(1,v[4],v[5],1)           # Params ckf(2)  ckf3(3)


  okuns <- dlmModReg(mod_data.tvp[,-1], addInt = TRUE, dV =dV, dW = dW, m0 = beta.vec)




}

#'Set parameter starting values

ckf4Guess <- -2.8913
ckf2guess <- -0.0495
ckf3guess <- 0.01942
ckf5guess <- -4.1757
ckf6guess <- -6.2466



params <- c(ckf4Guess,
            ckf5guess,
            ckf6guess,
            ckf2guess,
            ckf3guess)

tvp_mod.mle <- dlmMLE(mod_data.tvp[,"dur"] , parm = params, build = buildSS)

tvp_mod <- buildSS(tvp_mod.mle$par)

tvp_filter <-  dlmFilter(mod_data$dur,tvp_mod)

Приведенный выше код работает, но выходные данные не верны. Я не определяю штаты должным образом. У кого-нибудь есть опыт построения dlms с регрессом мутлвирата в R?

Answer 1

Я думаю, что нашел решение - мне удалось воссоздать оценки в статье, которая оценивает эту модель, используя Eviews (также проверил это, используя Eviews).

    #--------------------------------------------------------------------------------------------------------------------------
# tvp model full model - dur = alpha*dur(-1)+ beta(dgdp-potential) + gamma*wages
#--------------------------------------------------------------------------------------------------------------------------

# Construct DLM

OkunsDLMfm <- dlm(


  FF = matrix(c(1,1,1,1),ncol = 4, byrow = TRUE),

  V = matrix(1),

  GG = matrix(c(1,0,0,0,
                0,1,0,0,
                0,0,1,0,
                0,0,0,1), ncol = 4, byrow = TRUE),

  W =  matrix(c(1,0,0,0,
                0,1,0,0,
                0,0,1,0,
                0,0,0,1), ncol = 4, byrow = TRUE),

  JFF = matrix(c(1,2,3,0),ncol = 4, byrow = TRUE),

  X = cbind(mod_data$dur_lag1,mod_data$d2lgdp, mod_data$d2lrulc_lag2), # lagged dep var, dgdp, wages.

  m0 = c(0,0,0,0),

  C0 = matrix(c(1e+07,0,0,0,
                0,1e+07,0,0,
                0,0,1e+07,0,
                0,0,0,1e+07), ncol = 4, byrow = TRUE)

)


buildOkunsFM <- function(p){

  V(OkunsDLMfm)  <- exp(p[2])

  GG(OkunsDLMfm)[1,1]  <- 1

  GG(OkunsDLMfm)[2,2]  <- 1

  GG(OkunsDLMfm)[3,3]  <- 1 

  GG(OkunsDLMfm)[4,4]  <- 1

  W(OkunsDLMfm)[1,1] <- exp(p[3])

  W(OkunsDLMfm)[2,2] <- 0

  W(OkunsDLMfm)[3,3] <- 0

  W(OkunsDLMfm)[4,4] <- exp(p[4])

  m0(OkunsDLMfm) <- c(0,0,0,p[1]*4)

  C0(OkunsDLMfm)[1,1] <- 1

  C0(OkunsDLMfm)[4,4] <- 5


  return(OkunsDLMfm)

}



okuns.estfm <-  dlmMLE(y = mod_data$dur, parm = c(-0.049,-1.4,-6,-5), build = buildOkunsFM)


OkunsDLM1fm <- buildOkunsFM(okuns.estfm$par)

Изменяющийся во времени уровень, оценка потенциального выхода, получается путем деления 4 элемента вектора состояния на второй * на отрицательный 1.

Не уверен, что это лучший способ указать DLM, но результаты модели очень близки к полученным (в пределах 0,01) результатов от использования Eviews. При этом очень открыт для любых других спецификаций.