RE: [TMB users] Vector AR(1) random effects

[Kasper Kristensen]

Note that densities from from namespace 'density' return the negative log likelihood. So code should be:

    nll += VECSCALE(UNSTRUCTURED_CORR(rho),sigma)(error);

---

From: tmb-...@googlegroups.com [tmb-...@googlegroups.com] on behalf of Jarle Tufto [jarle...@gmail.com]
Sent: Monday, July 03, 2017 2:15 PM
To: TMB Users
Subject: [TMB users] Vector AR(1) random effects

Hi all -

I'm trying to implement random effects having a 2-dimensional vector AR(1) structure in TMB, that is,

  x_t = Phi x_{t-1} + w_t.

where w_t ~ N(0, Sigma).  This is different from 

  http://kaskr.github.io/adcomp/classdensity_1_1AR1__t.html 

which assumes that Phi is diagonal with equal elements along the diagonal.  The observations y are equal to x plus iid normal errors.

I want to fit the model under the assumption that the process is stationary.  I therefore parameterize the 2x2 matrix Phi in terms of the two eigenvalues lambda_1 and lambda_2 (or more precisely logit((lambda_i+1)/2) and the angles theta_1 and theta_2 of the two unit eigenvectors.  So this restricts Phi to have non-complex eigenvalues.  Sigma is parameterized in the usual way.  

The following code (also available at https://github.com/jtufto/tmb-var1.git) compiles and runs but the model doesn't converge as the inner optimization stops after 

233 iterations with the error 

  Error in iterate(par) : Newton dropout because inner gradient too steep.

I suspect that there is some error somewhere in my thinking or in my code as I would expect this be easily identifiable with a 2x500 observations of y.

Any comments will be greatly appreciated.

Jarle

#include <TMB.hpp>

template<class Type>

Type objective_function<Type>::operator() ()

{

  using namespace density;

  DATA_MATRIX(y);

  

  PARAMETER_MATRIX(x);

  PARAMETER(logsdy); Type sdy = exp(logsdy);

  

  // Sigma parameterized in terms of log of the two sd's and logit2 of rho

  PARAMETER_VECTOR(log_sigma); vector<Type> sigma = exp(log_sigma);

  PARAMETER_VECTOR(rho);

  

  // Phi parameterized in terms

  PARAMETER_VECTOR(theta);

  PARAMETER_VECTOR(logit_lambda); vector<Type> lambda = 2/(1+exp(-logit_lambda))-1; // eigenvalues

  

  matrix<Type> eigvec(2,2);

  eigvec(0,0) = cos(theta(0)); eigvec(0,1) = cos(theta(1));

  eigvec(1,0) = sin(theta(0)); eigvec(1,1) = sin(theta(1));

  

  matrix<Type> D(2,2);

  D(0,0) = lambda(0); D(0,1) = 0;

  D(1,0) = 0;         D(1,1) = lambda(1);

  

  matrix<Type> Phi = eigvec * D * eigvec.inverse();

  

  Type nll = 0;

  // Density of VAR(1) process random effects x

  for (int t=1; t<x.rows(); t++) {

    vector<Type> xt = x.row(t);

    vector<Type> xtminus1 = x.row(t-1);

    vector<Type> Phix = Phi*xtminus1;

    vector<Type> error = xt - Phix;

    nll -= VECSCALE(UNSTRUCTURED_CORR(rho),sigma)(error);

  }

  // Conditional density of observations y given x

  for (int t=0; t<y.rows(); t++) {

    for (int i=0; i<2; i++) {

      nll -= dnorm(y(t,i), x(t+x.rows()-y.rows(),i), sdy, true);

    }

  }

  

  return nll;

}

# setup Phi with real eigenvalues smaller than 1 in absolute value

theta <- c(0,pi/2)

lambda <- c(.5,.2)

vec <- matrix(c(cos(theta),sin(theta)),2,2,byrow=FALSE)

Phi <- vec %*% diag(lambda) %*% solve(vec)

Phi

eigen(Phi)

# Sigma

rho <- -.5

sigma <- c(2,3)

Rho <- diag(2)

Rho[2,1] <- Rho[1,2] <- rho

Sigma <- diag(sigma) %*% Rho %*% diag(sigma)

Sigma

# Stationary covariance matrix

Gamma0 <- solve(diag(4)-kronecker(Phi,Phi), as.vector(Sigma))

Gamma0 <- matrix(Gamma0,2,2)

Gamma0

# Simulate some data

n <- 500

x <- matrix(NA,n,2)

x[1,] <- mvtnorm::rmvnorm(1,sigma=Gamma0)

for (t in 2:n) 

  x[t,] <- Phi %*% x[t-1,] + t(mvtnorm::rmvnorm(1,sigma=Sigma))

# Check stationary covariance against theoretical value

var(x)

Gamma0

# Add some noise to get the data

sdy <- .1

y <- x + rnorm(2*n, sd=sdy)

#matplot(1:n, y, type="l",col=1)

# Fit the model

library(TMB)

compile("var1.cpp")

dyn.load(dynlib("var1"))

data <- list(y=y)

burnin <- 10

parameters <- list(

  x = matrix(0,nrow(y)+burnin, 2),

  logsdy = log(sdy),

  log_sigma = log(sigma),

  rho = rho,

  theta = c(0, pi/2), 

  logit_lambda = qlogis((lambda+1)/2)

)

obj <- MakeADFun(data,parameters,random="x")

obj$method="BFGS"

opt <- do.call(optim,obj)

sdreport(obj)

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/tmb-users/56c9ea1a-ca9a-497a-99c4-ecbbd33a4fbf%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Jarle Tufto]

Excellent! That certainly makes a difference!