/
simple.R
33 lines (31 loc) · 900 Bytes
/
simple.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
require(TMB)
dyn.load(dynlib("simple"))
## Test data
set.seed(123)
y <- rep(1900:2010,each=2)
year <- factor(y)
quarter <- factor(rep(1:4,length.out=length(year)))
period <- factor((y > mean(y))+1)
## Random year+quarter effect, fixed period effect:
B <- model.matrix(~year+quarter-1)
A <- model.matrix(~period-1)
B <- as(B,"dgTMatrix")
A <- as(A,"dgTMatrix")
u <- rnorm(ncol(B)) ## logsdu=0
beta <- rnorm(ncol(A))*100
eps <- rnorm(nrow(B),sd=1) ## logsd0=0
x <- as.numeric(A%*%beta+B%*%u+eps)
## Fit model
map <- NULL
obj <- MakeADFun(data=list(x=x,B=B,A=A),
parameters=list(u=u*0+.1,beta=beta*0+.1,logsdu=.1,logsd0=0.1),
random="u"
)
newtonOption(smartsearch=FALSE)
obj$fn(obj$par)
obj$gr(obj$par)
obj$control <- list(parscale=obj$par*0+1e-1,trace=10)
obj$hessian <- TRUE
opt <- do.call("optim",obj)
##summary(as.mle2(opt))
summary(as.mle(opt))