Seemingly Unrelated Regressions Time Series

[Sam Weiss]

Hi,

I'm trying to model a two time series as random walks with correlated error terms. However, the sampler doesn't complete its warmup. Any help would be appreciated. Thanks.

N is the number of data points and M is the number of series (in this case it is two). 

data {

int<lower=0> N;

int<lower=1> M;

vector[M] y[N];

}

parameters { 

vector[M] mu;

real<lower=0> sigma[2];

real<lower=-1,upper=1> rho;

}

transformed parameters {

cov_matrix[2] Sigma;

Sigma[1,2] <- rho* sigma[1] * sigma[2];

Sigma[2,1] <- Sigma[1,2];

Sigma[1,1] <- square(sigma[1]);

Sigma[2,2] <- square(sigma[2]);

}

model{

for(n in 2:N){

y[n]~multi_normal(y[n-1],Sigma);

}

}

[Bob Carpenter]

When you say "doesn't complete its warmup", what exactly happens?
I would expect the model below to converge pretty quickly. I don't
see anything wrong with the way it's written.

And what's the data? Have you tried simulating data according to
the model and fitting it?

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Sam Weiss]

Hi Bob,

Thanks for responding. I'm trying this on monthly car sales. 

R will get stuck with this current read out:

"

TRANSLATING MODEL 'dlm.monthly.dummies.2' FROM Stan CODE TO C++ CODE NOW.

COMPILING THE C++ CODE FOR MODEL 'dlm.monthly.dummies.2' NOW.

SAMPLING FOR MODEL 'dlm.monthly.dummies.2' NOW (CHAIN 1).

Iteration:   1 / 1000 [  0%]  (Warmup)

"

Enclosed is that data and below is the code I used. :

dlm.monthly.dummies.2="

data {

int<lower=0> N;

int<lower=1> M;

vector[M] y[N];

}

parameters { 

vector[M] mu;

real<lower=0> sigma[2];

real<lower=-1,upper=1> rho;

}

transformed parameters {

cov_matrix[2] Sigma;

Sigma[1,2] <- rho* sigma[1] * sigma[2];

Sigma[2,1] <- Sigma[1,2];

Sigma[1,1] <- square(sigma[1]);

Sigma[2,2] <- square(sigma[2]);

}

model{

for(n in 2:N){

y[n]~multi_normal(y[n-1],Sigma);

}

}

"

data=list(N=length(btn1[,1]),M=2, y=sales)

fit <- stan(model_code = dlm.monthly.dummies.2, data = data,  iter = 1000, chains = 4)

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ZkAaXcT539M/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Herra Huu]

Remove the unnecessary mu parameter vector from the model and it will work fine.

[Bob Carpenter]

Good catch.

The problem with the mu parameter with no distribution
specified is that Stan treats it as uniform on (-infinity, infinity),
which leads to an improper posterior. Stan requires proper
posteriors to sample. One way to guarantee this is with proper
priors, but Stan can handle improper priors as long as the posterior
is proper. There's a chapter discussing why such posteriors are
problematic in the manual.

- Bob

On Jun 13, 2014, at 4:53 PM, Herra Huu <heikki.uu...@gmail.com> wrote:

> Remove the unnecessary mu parameter vector from the model and it will work fine.
>

[Sam Weiss]

I see. thanks for the tip. It did the trick!

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ZkAaXcT539M/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Sam Weiss]

Hi again. I'm running the code below (same data) but its taking over 24 hours to complete... I'd like to do this with more than 2 time series as well. Is the code misspecified or inefficiently writen? 

Thanks,

library(rstan) 

set_cppo("fast") 

code1=

  "data {

int<lower=0> N;

int<lower=1> K;

vector[K] y[N];

}

parameters { 

corr_matrix[K] OmegaV;

vector<lower=0>[K] tauV;

corr_matrix[K] OmegaW;

vector<lower=0>[K] tauW;

vector[K] theta[N];

}

transformed parameters{

matrix[K,K] SigmaV;

matrix[K,K] SigmaW;

cov_matrix[2] Sigmavprior;

SigmaV <- diag_matrix(tauV) * OmegaV * diag_matrix(tauV);

SigmaW <- diag_matrix(tauW) * OmegaW * diag_matrix(tauW);

Sigmavprior[1,1] <- 100;

Sigmavprior[2,2]<- 100;

Sigmavprior[1,2] <- 0;

Sigmavprior[2,1] <- 0;  

}

model{

//prior

theta[1]~multi_normal(y[1], Sigmavprior);

//likelihood

for(n in 2:N){

theta[n]~multi_normal(theta[n-1],SigmaW);

y[n]~multi_normal(theta[n],SigmaV);

}

}

"

data=list(N=length(btn1[,1]),K=2, y=btn1[,c("American.Honda","Toyota")])

fit1 <- stan(model_code = code1, data = data,  iter = 1000, chains = 4)

[Bob Carpenter]

On Jun 23, 2014, at 7:15 PM, Sam Weiss <samc...@gmail.com> wrote:

> Hi again. I'm running the code below (same data) but its taking over 24 hours to complete... I'd like to do this with more than 2 time series as well. Is the code misspecified or inefficiently writen?
>
> Thanks,
>
>
> library(rstan)
> set_cppo("fast")
>
>
> code1=
> "data {
> int<lower=0> N;
> int<lower=1> K;
> vector[K] y[N];
>
> }
> parameters {
> corr_matrix[K] OmegaV;
> vector<lower=0>[K] tauV;
> corr_matrix[K] OmegaW;
> vector<lower=0>[K] tauW;
> vector[K] theta[N];
>
> }
> transformed parameters{
> matrix[K,K] SigmaV;
> matrix[K,K] SigmaW;
> cov_matrix[2] Sigmavprior;
>
>
> SigmaV <- diag_matrix(tauV) * OmegaV * diag_matrix(tauV);
>
> SigmaW <- diag_matrix(tauW) * OmegaW * diag_matrix(tauW);

There are more efficient scaling operations. diag_pre_multiply()
and diag_post_multiply() in Stan 2.2. In Stan 2.3 (out any
minute, really), there will also be a direct quadratic form that
will do this all at once.

>
>
> Sigmavprior[1,1] <- 100;
> Sigmavprior[2,2]<- 100;
> Sigmavprior[1,2] <- 0;
> Sigmavprior[2,1] <- 0;
>
> }

If you have anything constant like Sigmavprior, it's much
more efficient to define it as transformed data. All variables
should basically be declared in the earliest block that's feasible.

Also, no point making a diagonal covariance matrix. You can just
sample the components independently with a vectorized univariate
normal, replacing

> theta[1]~multi_normal(y[1], Sigmavprior);

with

theta[1] ~ multi_normal(y[1],100);

Now having said all this, you probably want to work directly
with Cholesky factors. You can do this one of two ways. Either
declare the Cholesky factors as variables (more efficient), or
Cholesky factorize a covariance matrix in the model itself and reuse
it with multi_normal_cholesky(). Alas, the Cholesky factor for
covariance matrix type is broken due to a misspelled error check in
2.2 (Ben --- this is why we try to be thorough in our signature
tests --- this one obviously didn't have enough testing); that's fixed
in 2.3 and there's also a Cholesky factor for correlation matrices type
and some explicit advice on how to use it in the manual.

> model{
> //prior
> theta[1]~multi_normal(y[1], Sigmavprior);
> //likelihood
> for(n in 2:N){
> theta[n]~multi_normal(theta[n-1],SigmaW);
> y[n]~multi_normal(theta[n],SigmaV);
>
> }
>

The implicit prior on covariances is improper, which can slow
down sampling --- are you seeing big excursions?

I like the first line written the other way around as
the equivalent

y[1] ~ multi_normal(theta[1], Sigmavprior);

to make it clear it's really part of the likelihood and to
make it clear you don't really have a prior here for theta[1].
Then you can add

theta[1] ~ normal(0,100);

if you want, but that prior's probably so weak as to have almost
no effect if your data are on a standardized scale.

- Bob

[Sam Weiss]

Hi bob,

Thanks again.

I got this error when I replaced Sigmavprior with:   theta[1] ~ multi_normal(y[1],100);

"

TRANSLATING MODEL 'code1' FROM Stan CODE TO C++ CODE NOW.

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  failed to parse Stan model 'code1' with error message:

EXPECTATION FAILURE LOCATION: file=input; line=27, column=34

theta[1]~multi_normal(y[1], 100);

                                 ^-- here

DIAGNOSTIC(S) FROM PARSER:

no matches for function name="multi_normal_log"

    arg 0 type=vector

    arg 1 type=vector

    arg 2 type=int

available function signatures for multi_normal_log:

0.  multi_normal_log(vector, vector, matrix) : real

unknown distribution=multi_normal

Parser expecting: <eps>

"

Will having a multi_normal(mu, 100) be the same as having a covariance matrix with 100 on the diagonal and zeros everywhere else?

Also, is this the correct forum to post these questions? 

[Daniel Lee]

Hi Sam,

Bob meant:

theta[1] ~ normal(y[1], 100);

instead of

theta[1] ~ multi_normal(y[1], 100);

On Mon, Jun 23, 2014 at 2:43 PM, Sam Weiss <samc...@gmail.com> wrote:

Hi bob,

Thanks again.

I got this error when I replaced Sigmavprior with:   theta[1] ~ multi_normal(y[1],100);

"

TRANSLATING MODEL 'code1' FROM Stan CODE TO C++ CODE NOW.

Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 

  failed to parse Stan model 'code1' with error message:

EXPECTATION FAILURE LOCATION: file=input; line=27, column=34

theta[1]~multi_normal(y[1], 100);

                                 ^-- here

DIAGNOSTIC(S) FROM PARSER:

no matches for function name="multi_normal_log"

    arg 0 type=vector

    arg 1 type=vector

    arg 2 type=int

available function signatures for multi_normal_log:

0.  multi_normal_log(vector, vector, matrix) : real

unknown distribution=multi_normal

Parser expecting: <eps>

"

Will having a multi_normal(mu, 100) be the same as having a covariance matrix with 100 on the diagonal and zeros everywhere else?

Also, is this the correct forum to post these questions? 

Yes. Please feel free to post questions to the Stan users list.

[Sam Weiss]

Hey Guys, 

I've successfully got my model to work! Very exciting.. Now I'd like to forecast my time series. I can't find any documentation to do this (no results come up when searching for "forecasting" in stan reference). Can someone suggest a place to look?

Thanks,

Sam

Also, code that worked is pasted below. It takes about 2 hours to complete...if anyone knows how to speed that up please let me know.

ode1=

  "data {

int<lower=0> N;

int<lower=1> K;

vector[K] y[N];

  int<lower=1>  S;

  int<lower=1>  dum[N]; // monthly dummies

}

parameters { 

corr_matrix[K] OmegaV;

vector<lower=0>[K] tauV;

corr_matrix[K] OmegaW;

vector<lower=0>[K] tauW;

vector[K] theta[N];

vector[K]  alpha[S];

vector[K] mu[N];

vector<lower=0>[K] sigma_mu;

}

transformed parameters{

matrix[K,K] L_betaV;

matrix[K,K] L_betaW;

L_betaV <- cholesky_decompose(OmegaV);

for (k1 in 1:K)

for (k2 in 1:k1)

L_betaV[k1,k2] <- tauV[k1] * L_betaV[k1,k2];

L_betaW <- cholesky_decompose(OmegaW);

for (k1 in 1:K)

for (k2 in 1:k1)

L_betaW[k1,k2] <- tauW[k1] * L_betaW[k1,k2];

}

model{

//prior

mu[1] ~ normal(0, 1.0e+5 );

theta[1]~normal(y[1], 100);

y[1]~multi_normal_cholesky(theta[1], L_betaV);

//likelihood

for(n in 2:N){

mu[n]~normal(mu[n-1], sigma_mu);

theta[n]~multi_normal_cholesky(mu[n]+theta[n-1],L_betaW);

y[n]~multi_normal_cholesky(alpha[ dum[n] ]  + theta[n],L_betaV);

}

}

"

[Bob Carpenter]

We talk about forecasting as "predictive inference". We should
add the "forecasting" terminology to aid in search. We talk about it in the
description of Bayesian inference in the intro and later on, and in a bunch of
examples, such as Gaussian processes and naive Bayes. I added a note
to the next manual issue to add an example for regression and time series:

https://github.com/stan-dev/stan/issues/668#issuecomment-47424210

There are two ways to do this.

Both involve adding new data variables for the number of items to
predict/forecast and any inputs necessary to do that prediction.
Then you can either:

1. Define new parameters for the predictions and give these parameters
predictions just like in the original model. This only works for continous
predictions.

2. Define new generated quantity variables for the predictions and
define them using random number generators in the generated quantities block.
This will work for discrete or continous predictions.

- Bob

[Sam Weiss]

This code doesn't work for a large number of variables (time series). When I run the code for two time series, it takes a long time to simulate (about an hour) and the results adequate. When I want to include something like 10 variables (in this case have 10 time series in my SUR), it takes a minute and gives values of all parameters around zero. 

Why does including more variables (in this case having a larger K) give inadequate results? (My thought is that the parameter space is too big / constrained and its not sampling properly.. Just a guess)

Thanks Again

Below is code: 

code1=

  "data {

int<lower=0> N1;

int<lower=0> N;

int<lower=0> N2;

int<lower=1> K;

vector[K] y1[N1];

int<lower=1>  S;

int<lower=1>  dum[N]; // monthly dummies

}

parameters { 

corr_matrix[K] OmegaV;

vector<lower=0>[K] tauV;

corr_matrix[K] OmegaW;

vector<lower=0>[K] tauW;

vector[K] theta[N];

vector[K]  alpha[S];

vector[K] mu[N];

vector<lower=0>[K] sigma_mu;

vector[K] y2[N2];

}

transformed parameters{

matrix[K,K] L_betaV;

matrix[K,K] L_betaW;

L_betaV <- cholesky_decompose(OmegaV);

for (k1 in 1:K)

for (k2 in 1:k1)

L_betaV[k1,k2] <- tauV[k1] * L_betaV[k1,k2];

L_betaW <- cholesky_decompose(OmegaW);

for (k1 in 1:K)

for (k2 in 1:k1)

L_betaW[k1,k2] <- tauW[k1] * L_betaW[k1,k2];

}

model{

vector[K] y[N];

for(n in 1:N1) y[n] <-y1[n];

for(n in 1:N2) y[N1+n] <-y2[n];

//prior

mu[1] ~ normal(y[2]-y[1], 10000);

theta[1]~normal(y[1], 5000);

y[1]~multi_normal_cholesky(theta[1], L_betaV);

//likelihood

for(n in 2:N){

mu[n]~normal(mu[n-1], sigma_mu);

theta[n]~multi_normal_cholesky(mu[n-1]+theta[n-1],L_betaW);

y[n]~multi_normal_cholesky(alpha[ dum[n] ]  + theta[n],L_betaV);

}

}

"

data1=list(N=length(data[,1])+12, N1=length(data[,1]),N2=12 ,K=10, y1=data.1[,c(1:10)], S=12, dum=rep(1:12,12)[1:(length(data[,1])+12)])

[Bob Carpenter]

Do you get warnings printed out as you go about failures
in the covariance matrix tests?

I'm hoping the problem is that our covariance matrix and
correlation matrix consistency checks are too stringent and
causing unnecessary failures. That's been the case before.

The underlying problem with more variables is that floating
point arithmetic just isn't that precise, and when you have
lots of operations (as in Cholesky factor multiplication
and solvers converting the constrained/unconstrained), the
results vary from what you'd get from infinite-precision
arithmetic in that we can lose symmetry or positive definiteness.

Another problem is if estimates get pushed up against boundaries
(here correlation = 1 or -1, sd = 0 or +infinity). This is
also very unstable arithmetically.

You might want to try again with Stan 2.4 --- we tried to clean
some of these things up, but we probably still have a ways to go.

- Bob

[Sam Weiss]

Hi Bob,

Tried it again with Stan 2.4. 

I do get one error at the beginning of each chain about the covariance matrix, it occurs only once. 

The distributions of some variables within a chain are all constant. 

Think changing the priors or including some could help?