Dopasowanie zmiennego w czasie współczynnika DLM

Chcę dopasować DLM ze zmiennymi w czasie współczynnikami, tj. Rozszerzeniem zwykłej regresji liniowej,

$y_t = \theta_1 + \theta_2x_2$.

Mam predyktor ($x_2$) i zmienną odpowiedzi ($y_t$), coroczne połowy ryb morskich i śródlądowych odpowiednio w latach 1950–2011. Chcę, aby postępował model regresji DLM,

$y_t = \theta_{t,1} + \theta_{t,2}x_t$

gdzie jest równanie ewolucji systemu

$\theta_t = G_t \theta_{t-1}$

od strony 43 Dynamicznych modeli liniowych z R autorstwa Petrisa i in.

Trochę kodowania tutaj,

fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4/fishdata.csv", header=T)
x <- fishdata$marinefao
    y <- fishdata$inlandfao

lmodel <- lm(y ~ x)
summary(lmodel)
plot(x, y)
abline(lmodel)

W tym przypadku bardziej odpowiednie są zmienne w czasie współczynniki modelu regresji. Śledzę jego przykład ze stron 121–125 i chcę zastosować to do moich danych. To jest kodowanie z przykładu

############ PAGE 123
require(dlm)

capm <- read.table("http://shazam.econ.ubc.ca/intro/P.txt", header=T)
capm.ts <- ts(capm, start = c(1978, 1), frequency = 12)
colnames(capm)
plot(capm.ts)
IBM <- capm.ts[, "IBM"]  - capm.ts[, "RKFREE"]
x <- capm.ts[, "MARKET"] - capm.ts[, "RKFREE"]
x
plot(x)
outLM <- lm(IBM ~ x)
outLM$coef
    acf(outLM$res)
qqnorm(outLM$res)
    sig <- var(outLM$res)
sig

mod <- dlmModReg(x,dV = sig, m0 = c(0, 1.5), C0 = diag(c(1e+07, 1)))
outF <- dlmFilter(IBM, mod)
outF$m
    plot(outF$m)
outF$m[ 1 + length(IBM), ]

########## PAGES 124-125
buildCapm <- function(u){
  dlmModReg(x, dV = exp(u[1]), dW = exp(u[2:3]))
}

outMLE <- dlmMLE(IBM, parm = rep(0,3), buildCapm)
exp(outMLE$par)
    outMLE
    outMLE$value
mod <- buildCapm(outMLE$par)
    outS <- dlmSmooth(IBM, mod)
    plot(dropFirst(outS$s))
outS$s

Chcę mieć możliwość kreślenia oszacowań wygładzania plot(dropFirst(outS$s))dla moich własnych danych, które mam problem z wykonaniem.

AKTUALIZACJA

Mogę teraz produkować te wykresy, ale nie sądzę, aby były prawidłowe.

fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4/fishdata.csv", header=T)
x <- as.numeric(fishdata$marinefao)
    y <- as.numeric(fishdata$inlandfao)
xts <- ts(x, start=c(1950,1), frequency=1)
xts
yts <- ts(y, start=c(1950,1), frequency=1)
yts

lmodel <- lm(yts ~ xts)
#################################################
require(dlm)
    buildCapm <- function(u){
  dlmModReg(xts, dV = exp(u[1]), dW = exp(u[2:3]))
}

outMLE <- dlmMLE(yts, parm = rep(0,3), buildCapm)
exp(outMLE$par)
        outMLE$value
mod <- buildCapm(outMLE$par)
        outS <- dlmSmooth(yts, mod)
        plot(dropFirst(outS$s))

> summary(outS$s); lmodel$coef
       V1              V2       
 Min.   :87.67   Min.   :1.445  
 1st Qu.:87.67   1st Qu.:1.924  
 Median :87.67   Median :3.803  
 Mean   :87.67   Mean   :4.084  
 3rd Qu.:87.67   3rd Qu.:6.244  
 Max.   :87.67   Max.   :7.853  
 (Intercept)          xts 
273858.30308      1.22505 

Oszacowanie wygładzania przecięcia (V1) jest dalekie od współczynnika regresji lm. Zakładam, że powinni być bliżej siebie.

Jaki jest dokładnie twój problem?

Jedyną pułapką, jaką znalazłem, jest to, że najwyraźniej

fishdata <- read.csv("http://dl.dropbox.com/s/4w0utkqdhqribl4,
                     fishdata.csv", header=T)

odczytuje dane jako liczby całkowite. Musiałem przekonwertować je na pływak,

x <- as.numeric(fishdata$marinefao)
y <- as.numeric(fishdata$inlandfao)

zanim mogłem wywołać funkcje dlm *.