9
$\begingroup$

I want to fit a DLM with time-varying coefficients, i.e. an extension to the usual linear regression,

$y_t = \theta_1 + \theta_2x_2$.

I have a predictor ($x_2$) and a response variable ($y_t$), marine & inland annual fish catches respectively from 1950 - 2011. I want the DLM regression model to follow,

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

where the system evolution equation is

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

from page 43 of Dynamic Linear Models With R by Petris et al.

Some coding here,

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)

Clearly time-varying coefficients of the regression model are more appropriate here. I follow his example from pages 121 - 125 and want to apply this to my own data. This is the coding from the example

############ 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

I want to be able to plot the smoothing estimates plot(dropFirst(outS$s)) for my own data, which I'm having trouble executing.

UPDATE

I can now produce these plots but I don't think they are correct.

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 

The intercept smoothing estimate (V1) is far from the lm regression coefficient. I assume they should be nearer to each other.

$\endgroup$
2
$\begingroup$

What is exactly your problem?

The only pitfall I found is that, apparently,

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

reads data as integers. I had to convert them to float,

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

before I could invoke the dlm* functions.

$\endgroup$
  • $\begingroup$ Thanks for your suggestions @F. Tusell; I have updated my question. The smoothing estimates produced are not close to the lmodel$coef estimates. I assume the plots are incorrect but I could be wrong. $\endgroup$ – hgeop Feb 9 '14 at 18:11
  • 1
    $\begingroup$ There is no reason to expect the smoothed estimates of slope and intercept to be close to the fixed betas in the linear regression. In particular, the slope should fluctuate wildly. $\endgroup$ – F. Tusell Feb 10 '14 at 8:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.