koyck model in stan

[Mark Leeds]

Hi to everyone on the STAN users mailng list : I'll try to be brief but I wanted to explain the history of my problem first. The model I'm trying to fit is called the koyck disributed lag and has  a long classical history. In fact, it was developed in 1951
in the Ph.D thesis of Koyck. Phoebus dhrymes actually wrote a whole book on it and there have been many, many papers on it because estimating it, even classically, was not easy before computers got powerful. This is because there
are similarities between the koyck model and the GLS model ( serial correlation ) that are not necessary to go into here.

But I'm not so sure about it's history as far as bayesian attempts to estimate it. The paper that attempts to estimate it using gibbs sampling is attached but I found an error in it ( equation 44 ) and the authors don't provide the explicit GIbbs sampling
 code that estimated the model so I'm not so sure about the correctness of the results.  Also, there is definitely an issue at the boundary when rho is near 1 so that could be the reason why I have had so many problems estimating it. The beginning of attached  paper provides a discussion of the gelman-meng example and how one can get ridges in the likelihood or a bi-modal likelihood depending on the true parameter values.  The actual koyck model is discussed on pages 19-22. ( although the discussion between pages 12-18 is also quite relevant ) and pages 23-24 talked about problems at the boundary and how to deal with them.

My history with the model is that I worked on the bayesian estimation of it about 3 years ago. I made three different attempts:

1) I tried Rjags. ( I have that code if anyone is interested ).

2) I wrote my own code using gibbs sampling. When that didn't work, I wrote a metropolis hastings algoriithm. That wasn't successful either. ( I have that R code also ).

3) I tried to use the R admit package. That package generates a mixture of student t-densities as the proposal for an MH algorithm. It was actually the most successful of the three attempts but when the rho parameters was near the boundaty,
it also had problems.

The main problem is that the model I'm estimating with real data would be expected to have estimates of rho in the (0.9 - 1.0) range so having problems near the boundary is a real problem.

So, I heard about STAN a month or so ago and, after reading about HMC and NUTS, I finally tried it tonight. I didn't get good parameter estimates but of course I could be doing something incorrectly. I've only just started reading the Jags manual but this model is actually not that different from an AR(1) because a koyck distributed lag reduces to something pretty similar to an AR(1). So, the AR(1) model in the STAN manual really helped a lot. 

From reading the manual and Ben's great postings, it seems that RSTAN assumes flat improper priors for parameters that aren't explicity given a prior. So, I think the model matches what is in the paper as far as the joint prior.

Also,  I realize that it's a lot of work to look at the paper and the model etc but, if  anyone is interested, I'd really appreciate it and maybe something interesting could come out of it. I originally emailed Andrew a month or so ago about making a
contest on estimating this model on his blog but then he mentioned STAN which I had zero knowledge of.  Thanks Andrew.

As far as how my 3 years ago journey ended,  I ended up estimating the model using optimization methods but I'd love to know if it actually can be estimated using any kind of MCMC approach. Plus, when you use classical methods, you
have to decide on what time window to use to re-update the model estimates etc. A bayesian approach would avoid the need for that.

Or maybe it's a rare example of a model that just cannot be estimated using an MCMC approach. Thank you to anyone who wants to take a look. If anyone is interested, then my suggestion would be to read pages 1-22 of the attached paper. It's well written and really clear. It's only the Gibbs sampling results that I'm skeptical of. Thank you very much and if anyone has any questions about the model, feel free to ask me. Also, one more piece of info: 3 years back, I did get help from Ed Greenberg as far as asking him about my GIbbs sampling R code  and he wasn't able to figure out why my Gibbs sampling R code wasn't working so  I'm pretty certain that it's not possible to estimate it using regular GIbbs sampling.

Finally, my thanks to all of the people involved in Stan and Rstan. The instructions for getting it to run are so well written that it was a breeze. I always have bad times with software installations and, aside from my dopey locale problem, it was a piece of cake getting STAN and Rstan to work on my linux machine. Also, although I only started reading it today, the STAN manual is incredible with all the examples and details and explanations. and the users list is amazingly educational also. Thank you all  very much.


                                                                                                                                                                                                                                                                              Mark

As an aside, if anyone knows anything that would help me understand NUTS more, please let me know. The authors, Andrew and Matt, are both well known and extremely clear in their writings but the connection between binary trees and their depth and using that information to decide when to stop the trajectory generation is something I'm just not getting.  Thanks again.

[Carlos]

Mmmm, why don't you assign priors? I have just changed your code as follows,

data {

    int<lower=0> N; // number of data points 

    real y[N];

    real x[N];   

}

parameters {

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

    real <lower=0> sigma;

    real beta;

}

transformed parameters{

    real <lower=0> sigmas;    

    sigmas <- sqrt(sigma);

}

model {

    for ( n in 2:N)

        y[n] ~ normal(rho * y[n-1] + beta * (1-rho) * x[n], sigmas);

    // Priors

    rho ~ uniform(0,1);

    beta ~ normal(0,100);

    sigma ~ inv_gamma(2,1);

}

and it works good... :) Without that, I got totally non-sense...

[Mark Leeds]

thank carlos. I'll try that but I was trying to mimic the priors in the paper which  are flat improper priors. I think the manual says that if you don't put any priors, then flat improper ones are assumed which is what I wanted to use.

but that's great that it works because I don't really care that much about what priors are used.  I'll play around with it and see how it works near the rho boundary. I really appreciate you taking the time to look at it. thanks again.

On Saturday, April 6, 2013 11:36:39 AM UTC-4, Carlos wrote:

Mmmm, why don't you assign priors? I have just changed your code as follows,

data {

    int<lower=0> N; // number of data points 

    real y[N];

    real x[N];   

}

parameters {

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

    real <lower=0> sigma;

    real beta;

}

transformed parameters{

    real <lower=0> sigmas;    

    sigmas <- sqrt(sigma);

}

model {

    for ( n in 2:N)

        y[n] ~ normal(rho * y[n-1] + beta * (1-rho) * x[n], sigmas);

    // Priors

    rho ~ uniform(0,1);

    beta ~ normal(0,100);

    sigma ~ inv_gamma(2,1);

}

and it works good... :) Without that, I got totally non-sense...


On Saturday, April 6, 2013 2:24:16 AM UTC-4, Mark Leeds wrote:

[Mark Leeds]

Hi Carlos: not that it matters, but one other thing: your uniform prior on rho is redundant given how STAN treats constraints on parameters. Ben explains all of this quite thoroughly in his posts at the top of the  STAN users mailing list.  I suggest checking Ben's posts out because they were really useful to me and not reading them could bite you when you when building larger, more complex models.

Thanks aagain and I'll get back to you about how well your code works.

[Ben Goodrich]

On Saturday, April 6, 2013 2:24:16 AM UTC-4, Mark Leeds wrote:

From reading the manual and Ben's great postings, it seems that RSTAN assumes flat improper priors for parameters that aren't explicity given a prior. So, I think the model matches what is in the paper as far as the joint prior.

Minor correction: The priors are uniform over their declared support, unless you bend the prior explicitly with a sampling statement in the model {} block such as beta ~ normal(0,1). But a uniform prior is proper if the support is bounded on both sides, which is the case here with rho in (0,1).

As to your model, the wall time would improve if you pulled the sampling statement out of the loop:

data {
  int<lower=0> N; // number of data points 
  real y[N];
  real x[N];
}

model {
 real mu[N];
 mu[0] <- y[0];

 for ( n in 2:N)

   mu[n] ~ rho * y[n-1] + beta * (1-rho) * x[n];

 y ~ normal(mu, sigma);
}

If you are quite sure that in your real example that rho is large, then you can change the lower bound of the support to .8 or so, and it will sample a bit better because the unconstrained parameter underlying rho will be closer to symmetric. Also, like Carlos said, some reasonable but weak priors would help too.

Ben

[Mark Leeds]

Hi Ben: Thanks for your comment. Just to make sure I understand: In the paper, the joint prior is f(beta, rho, sigma) ~ 1/sigma ( I'm not sure what is assumed individually to get that joint )  so if I do

A) lower=0, upper=1 for rho

and

B) nothing for the others,

am i mimicing the paper exactly ? I haven't had time to play around yet but I'd like to be exact, atleast initially, and see if I can get reasonable estimates for rho reasonably near 1.0.

Also, I'd rather not put a normal prior on beta because I basically have zero information for beta so I don't want to impose information. Same for sigma. rho has to be between 0 and 1 for stationarity ( actually -1 and + 1, but i know
it's positive )  so a uniform prior between 0 and 1 is fine. But , when you said (0.8,1) for rho,  that would also be true for say 0.7, 0.6 etc also, right ? Thanks a lot.

[Ben Goodrich]

On Saturday, April 6, 2013 9:35:00 PM UTC-4, Mark Leeds wrote:

Hi Ben: Thanks for your comment. Just to make sure I understand: In the paper, the joint prior is f(beta, rho, sigma) ~ 1/sigma ( I'm not sure what is assumed individually to get that joint ) 

I don't know what assumptions yield that joint either. My guess is such a prior puts mass on rho outside of (-1,1) even though the model would be non-stationary.
 

so if I do

A) lower=0, upper=1 for rho

and

B) nothing for the others,

am i mimicing the paper exactly ?

No, even though by "nothing" you are still doing something when you restrict sigma to be positive. Possibly if you define rho as a transformed parameter, such as tanh(theta) to be on (-1,1) or logit(theta) to be on (0,1) you could construct an improper prior for f(beta, theta, sigma).
 

I haven't had time to play around yet but I'd like to be exact, atleast initially, and see if I can get reasonable estimates for rho reasonably near 1.0.

I think the "near 1.0" problem is manageable. For instance, you could define rho as the transformed parameter theta^(1/K) for some theta in (0,1) and sufficiently large K, in which case theta would be sufficiently below 1.0.
 

Also, I'd rather not put a normal prior on beta because I basically have zero information for beta so I don't want to impose information. Same for sigma. rho has to be between 0 and 1 for stationarity ( actually -1 and + 1, but i know it's positive )

Sometimes there are situations where MCMC should work in principle without prior information but does not work in practice for particular starting values, seed, etc. So, if you just put some weak priors on things, it can work better.
 

  so a uniform prior between 0 and 1 is fine. But , when you said (0.8,1) for rho,  that would also be true for say 0.7, 0.6 etc also, right ? Thanks a lot.

Yeah, whatever lower bound you want is fine.

Another approach would go like this

parameters {

  real <lower=0, upper=1> rho;
  real <lower = 0> sigma;

  
  real beta1mrho; // equals beta * (1 - rho)

}
model {

 real mu[N];
 mu[0] <- y[0];

 for ( n in 2:N)
   mu[n] ~ rho * y[n-1] + beta * (1-rho) * x[n];

 y ~ normal(mu, sigma);

}

generated quantities {
  real beta;
  beta <- beta1mrho / (1 - rho)
}

That might get you better geometry.

Ben

[Mark Leeds]

Thanks Ben for all of your great help. Hopefully I'll get some time to play around tomorrow ( sunday ) and see how things go. Of course, I'll report back and let you  know the results. and yes, the joint prior in the paper is quite
strange to me also.

[Ben Goodrich]

Sounds great. In my previous post I meant to write

mu[n] ~ rho * y[n-1] + beta1mrho * x[n];

rather than

mu[n] ~ rho * y[n-1] + beta * (1-rho) * x[n];

Ben

[Mark Leeds]

thanks ben. I was confused a little by the original code but was afraid it was my lack of experience in STAN so I didn't ask. I'll get back to you but not exactly sure when. hopefully sooner than later.

also, you are probably already aware of this but there's a lot of ADL work in political science ( nathaniel beck is one name that comes to my mind ).  the koyck model is kind of a special case of an ADL  so you may be interested in checking it out that part of the poly sci literature. I have many papers so, if you're interested, I can dig them up and send you some titles. I'm not a political science person but I found that some of the political science literature gave beautiful expositions of what's really going on with the koyck distributed lag. nice to meet you.

[Mark Leeds]

Hi Carlos: I'm sorry to bother you but I just got back into this since we spoke and I can't find a set of parameters that works well even with the priors you included in your last email. I'm not claiming
that priors can't help and I haven't changed any of the code to benefit from Ben's comments ( all I did was add your priors ) but the attached code doesn't provide good estimates for the set of parameters
given. At your absolute convenience,  could you put your parameters in the code attached and run it and let me know what you get.  Thanks a lot.

On Saturday, April 6, 2013 11:36:39 AM UTC-4, Carlos wrote:

Mmmm, why don't you assign priors? I have just changed your code as follows,

data {

    int<lower=0> N; // number of data points 

    real y[N];

    real x[N];   

}

parameters {

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

    real <lower=0> sigma;

    real beta;

}

transformed parameters{

    real <lower=0> sigmas;    

    sigmas <- sqrt(sigma);

}

model {

    for ( n in 2:N)

        y[n] ~ normal(rho * y[n-1] + beta * (1-rho) * x[n], sigmas);

    // Priors

    rho ~ uniform(0,1);

    beta ~ normal(0,100);

    sigma ~ inv_gamma(2,1);

}

and it works good... :) Without that, I got totally non-sense...

On Saturday, April 6, 2013 2:24:16 AM UTC-4, Mark Leeds wrote:

[Bob Carpenter]

I just wrote a first-order Koyck model in Stan along with
a simulator in R. It all seems to work without a hitch.
Here are the links to the code:

R Simulator:
https://github.com/stan-dev/stan/blob/master/src/models/misc/garch/koyck-sim.R

Stan Model:
https://github.com/stan-dev/stan/blob/master/src/models/misc/garch/koyck.stan

The model (which is very much like Mark Leeds's) is short enough to dump
out in the mail:

// geometric lag time-series (Koyck 1951)
//

data {
int<lower=0> T; // number of time points
real y[T]; // output at time t
real x[T]; // predictor for time t
}
parameters {
real alpha; // intercept
real beta; // slope
real <lower=0, upper=1> lambda; // lag
real <lower=0> sigma; // noise scale
}
model {
alpha ~ cauchy(0,5);
beta ~ cauchy(0,5);
lambda ~ uniform(0,1);
sigma ~ cauchy(0,5);
for (t in 2:T)
y[t] ~ normal(alpha + beta * x[t] + lambda * y[t-1],
sigma);
}

I followed the parameterization in:

http://en.wikipedia.org/wiki/Distributed_lag

I used proper priors for everything, but fairly diffuse ones.
You could eliminate the uniform one, but having it isn't going
to slow things down noticeably.

One difference from Mike Leeds's Stan model is that I used an
intercept term alpha. Without an intercept, the outputs should
be centered (this could be done in Stan's transformed data block).

Here's the simulator in R:

library("rstan");

T <- 200;
x <- rnorm(T,0,1);
alpha <- -1.25;
beta <- 0.75;
lambda <- 0.6;
sigma <- 0.5;

y <- rep(NA,T);
y[1] <- rnorm(1,alpha + beta * x[1],sigma);
for (t in 2:T)
y[t] <- rnorm(1,alpha + beta * x[t] + lambda * y[t-1], sigma);

fit <- stan(file="koyck.stan", data=c("T","y","x"), iter=500, chains=4, fit=fit);


And here's the output, which is pretty much instant after the mode
compiles:

Inference for Stan model: koyck.
4 chains: each with iter=500; warmup=250; thin=1; 500 iterations saved.

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha -1.18 0.01 0.11 -1.38 -1.24 -1.18 -1.10 -0.97 287 1.01
beta 0.73 0.00 0.03 0.66 0.70 0.73 0.75 0.79 312 1.01
lambda 0.63 0.00 0.03 0.57 0.61 0.63 0.65 0.69 275 1.01
sigma 0.49 0.00 0.02 0.45 0.47 0.49 0.51 0.54 277 1.01
lp__ 39.73 0.09 1.36 36.17 39.01 39.99 40.78 41.51 256 1.00


I also tried a lambda=0.98 setting with 10K data points and
everything else the same, and here's the fit:

mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
alpha -1.269 0.003 0.058 -1.386 -1.312 -1.270 -1.229 -1.156 429 1.004
beta 0.750 0.000 0.005 0.741 0.747 0.750 0.753 0.759 1183 1.002
lambda 0.980 0.000 0.001 0.978 0.979 0.980 0.980 0.981 428 1.005
sigma 0.504 0.000 0.003 0.497 0.501 0.504 0.506 0.511 1065 1.003
lp__ 1857.899 0.063 1.353 1854.309 1857.265 1858.221 1858.858 1859.628 465 1.011

So I think it's safe to say that Stan can fit the Koyck model :-)

It's slower with a very high lambda, but NUTS's adaptation still
figures it out. The model could be made quite a bit more efficient
with vectorization.

- Bob

On 4/8/13 12:07 PM, Mark Leeds wrote:
> Hi Carlos: I'm sorry to bother you but I just got back into this since we spoke and I can't find a set of parameters
> that works well even with the priors you included in your last email. I'm not claiming
> that priors can't help and I haven't changed any of the code to benefit from Ben's comments

...

[Ben Goodrich]

Bob just created an example with one of these models that works with weakly informative priors

library("rstan");
T <- 200;
x <- rnorm(T,0,1);
alpha <- -1.25;
beta <- 0.75;
lambda <- 0.6;
sigma <- 0.5;

y <- rep(NA,T);
y[1] <- rnorm(1,alpha + beta * x[1],sigma);
for (t in 2:T)
  y[t] <- rnorm(1,alpha + beta * x[t] + lambda * y[t-1], sigma);

code <-

"

data {
  int<lower=0> T;   // number of time points
  real y[T];        // output at time t
  real x[T];        // predictor for time t
}
parameters {
  real alpha;                       // intercept
  real beta;                        // slope
  real <lower=0, upper=1> lambda;   // lag
  real <lower=0> sigma;             // noise scale
}
model {
  alpha ~ cauchy(0,5);
  beta ~ cauchy(0,5);

  // implicit: lambda ~ uniform(0,1);

  sigma ~ cauchy(0,5);
  for (t in 2:T)
    y[t] ~ normal(alpha +  beta * x[t] +  lambda * y[t-1],
                  sigma);
}
"

fit <- stan(model_code = code, data=c("T","y","x"), iter=500, chains=4);
fit

Ben

[Mark Leeds]

wow bob. thanks. I'll have to modify it for my case but that's great. I saw an example I think in jackman where he does centering. is that talked about anywhere in the literature ? I don't know if that's
ultimately what my problem was when I used rjags but I'll

b) take your code and  make sure I understand it

b) make a new koyck model with no alpha and centering that hopefully works.

c) put it somewhere on github and send out the links.

d) then I'll see if that same code also fixes the problems I was having in rjags. I'm still curious if the rjags problem is the jags algorithm or the centering issue ?

thanks so much. it was very nice of you to take the time to write the full model !!!!!

On Saturday, April 6, 2013 2:24:16 AM UTC-4, Mark Leeds wrote:

[Mark Leeds]

pages 270-272 in jackman's text "bayesian analysis for the social sciences" explain why centering in MCMC can be helpful. if anyone doesn't have the text and wants those pages, just contact
me off list and I can scan them and send them.  markl...@gmail.com.

On Saturday, April 6, 2013 2:24:16 AM UTC-4, Mark Leeds wrote:

[Mark Leeds]

Here are some initial comments/results concerning Bob's koyck model:

1) Unfortunately, the koyck model in the wiki is incorrect. They just wrote down the reduced form but in that reduced form, there should be a (1-lambda) multiplying beta. There are
two different versions  of the reduced form of the koyck distributed model but neither of them are the one in the wiki. The two forms result from the particular assumption put on the
initial error term in the long form. I generally refer to them ( in my lingo ) as the depooter model and the baltagi model ( without intercepts )

1) depooter model ( one from earlier attachment on page 19 of adl_gibbs ):   y_t = lambda* y_t-1 + beta(1-lambda)*X_t + epsilon_t ( replacing rho in the paper with lambda )

2) baltagi ( 6.10 on page 136 of attached ) : y_t = lambda*y_t-1 + beta*(1-lambda)*X_t + epsilon_t - lambda*epsilon_t-1   

I tend to work with depooter's  but I've also estimated 2) using numerical optimization in R.

Even though the wiki model doesn't match 1) or 2), I played with Bob's model in the following ways:

A) first I took alpha out and this had no negative effects. The estimates remained fine regardless of the value of lambda.

B) then I went ahead and took the model in A) and added the (1-lambda) term to the X_t term in order to make it the depooter model. This is where problems occur.  When
lambda is in the higher range, 0.7-0.95, the estimates of beta get worse as lambda gets closer to 1.0. This is because, the coefficient (1-lambda)*beta becomes less
indentifiable as lambda gets closer to 1.0 ( this is what is discussed in adl_gibbs.pdf on pages 23-24 ).  So, atleast now I know that the problem is due to the fact that lambda is involved
in both coefficients. I haven't tried centering yet  and that will be next effort but I have one sort of related question before I work on that.

In the classical formulation, ( where I estimate the model numerical using gradients or linear regression ), then if I estimated the Depooter model using Bob's original formulation
and then just divided the model beta estimate, namely beta hat,  by (1-lambda hat) to get a betastar hat, that's a legitimate "trick" to avoid dealing with the parameter being involved in both
coefficients. Bit this is only a valid "trick" when estimating it classically I think.

In the bayesian setting, that's cheating right, because the original model isn't being told that there's a (1-lambda) term on the beta being estimated in the original model ?
I just  want to make sure that doing that is illegal statistically in the bayesian formulation.  If it's legal, then I can use Bob's formulation to back out the value of beta hat
in my formulation. But that would seem illegal to me because one is making an assumption about the coeffs of the model and not telling  the model about that assumption during estimation ?
Anyway, I'm going to try that next and see what happens: Estimate the model using Bob's formulation but back out the estimate implied by my formulation and see if it's correct ?
 
Assuming that I'm correct about that being cheating, my next step after that will be centering. The Stan code and R code for the Depooter model without alpha is provided below.
I tried different "triplets of estimates" , using the same near 1.0 value of lambda each time and lowering it away from 1.0 using the next triplet to see when good estimates of beta  start to occur.
Thanks for any insights regarding my cheating question.

R CODE
#==========================================================================================================

library("rstan");

T <- 200;
x <- rnorm(T,0,1);

#=========================================
# NOT OKAY
#=========================================
#beta <- 1.0;
#lambda <- 0.95
#sigma <- 1.0

#beta <- 2.6;
#lambda <- 0.95
#sigma <- 1.0

#beta <- 0.5;
#lambda <- 0.95
#sigma <- 1.0
#====================================
# NOT OKAY
#====================================
#beta <- 1.0;
#lambda <- 0.90
#sigma <- 1.0

#beta <- 2.6;
#lambda <- 0.90
#sigma <- 1.0

#beta <- 0.5;
#lambda <- 0.90
#sigma <- 1.0
#========================================
#NOT OKAY
#=========================================
#beta <- 1.0;
#lambda <- 0.85
#sigma <- 1.0

#beta <- 2.6;
#lambda <- 0.85
#sigma <- 1.0

#beta <- 0.5;
#lambda <- 0.85
#sigma <- 1.0
#====================================
# NOT OKAY
#====================================
#beta <- 1.0;
#lambda <- 0.80
#sigma <- 1.0

#beta <- 2.6;
#lambda <- 0.80
#sigma <- 1.0

#beta <- 0.5;
#lambda <- 0.80
#sigma <- 1.0
#====================================
# NOT OKAY
#====================================
#beta <- 1.0;
#lambda <- 0.75
#sigma <- 1.0

#beta <- 2.6;
#lambda <- 0.75
#sigma <- 1.0

#beta <- 0.5;
#lambda <- 0.75
#sigma <- 1.0
#====================================
# ALMOST OKAY
#====================================
# OKAY
#beta <- 1.0;
#lambda <- 0.70
#sigma <- 1.0

# OKAY
#beta <- 2.6;
#lambda <- 0.70
#sigma <- 1.0

# STILL NOT OKAY
#beta <- 0.5;
#lambda <- 0.70
#sigma <- 1.0
#=============================================
# BETTER BUT NOT QUITE
#beta <- 0.5;
#lambda <- 0.65
#sigma <- 1.0
#=============================================
# OKAY !!!!!!
beta <- 0.5;
lambda <- 0.60
sigma <- 1.0
#===========================================
y <- rep(NA,T);

y[1] <- rnorm(1,beta * (1 - lambda) * x[1],sigma);

for (t in 2:T)

  y[t] <- rnorm(1, beta * (1 - lambda) * x[t] + lambda * y[t-1], sigma);

fit <- stan(file="koyck_nalpha_oneminus.stan", data=c("T","y","x"), iter=1000, chains=4);
print(fit)

#===========================================================================================================
# STAN CODE
#===========================================================================================================

data {
  int<lower=0> T; // number of time points
  real y[T]; // output at time t
  real x[T]; // predictor for time t
}
parameters {

  real beta; // slope
  real <lower=0, upper=1> lambda; // lag
  real <lower=0> sigma; // noise scale
}
model {

  beta ~ cauchy(0,5);

  lambda ~ uniform(0,1);
  sigma ~ cauchy(0,5);
  for (t in 2:T)

    y[t] ~ normal((1-lambda) * beta * x[t] + lambda * y[t-1],sigma);

[Andrew Gelman]

All I can say is, once this is done, Mark and Bob should publish a paper on this.  It will be helpful to the applied field, also good for Stan to start getting some publications.

A

--
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/groups/opt_out.
 
 
<adl_baltagi.pdf>

[Bob Carpenter]

Many of these time-series models assume that the
outputs are centered so that there's no intercept term.

Just think of basic regression. If you have

y[n] = alpha + beta * x[n] + epsilon[n]

and

epsilon[n] ~ normal(0,sigma)

then if alpha is estimated as non-zero, then the slope-only
model is going to be less effective,

y[n] = beta2 * x[n] + epsilon2[n]

and the slope beta2 won't match the original slope beta
because it needs to try to take up some of the work
formerly done by alpha. And the scale of the error
distribution for epsilon2 will be larger.

Now if you take y[n] and subtract the average value
(which is going to be alpha = mean(y) in a basic linear regression),

y2[n] <- y[n] - mean(y);

then the regression

y2[n] = beta * x[n] + epsilon[n];

I don't know what a good reference is for centering. It should
be discussed heavily within all these time series models because
they do it all the time.

There's a bit of discussion of completely standardizing
(y -> (y - mean(y))/sd(y)), which both centers and rescales.

- Bob

[Andrew Gelman]

We discuss centering and standardizing (dividing by 2 sd's, actually) in chapter 3 or 4 (I don't remember which) of ARM.

[Bob Carpenter]

On 4/8/13 7:41 PM, Andrew Gelman wrote:
> We discuss centering and standardizing (dividing by 2 sd's, actually) in chapter 3 or 4 (I don't remember which) of ARM.

By which Andrew means his applied regression modeling (hence the ARM)
book with Jennifer Hill.

- Bob

[Ben Goodrich]

On Monday, April 8, 2013 7:10:32 PM UTC-4, Mark Leeds wrote:

B) then I went ahead and took the model in A) and added the (1-lambda) term to the X_t term in order to make it the depooter model. This is where problems occur.  When
lambda is in the higher range, 0.7-0.95, the estimates of beta get worse as lambda gets closer to 1.0. This is because, the coefficient (1-lambda)*beta becomes less
indentifiable as lambda gets closer to 1.0 ( this is what is discussed in adl_gibbs.pdf on pages 23-24 ).  So, atleast now I know that the problem is due to the fact that lambda is involved
in both coefficients. I haven't tried centering yet  and that will be next effort but I have one sort of related question before I work on that.

In the classical formulation, ( where I estimate the model numerical using gradients or linear regression ), then if I estimated the Depooter model using Bob's original formulation
and then just divided the model beta estimate, namely beta hat,  by (1-lambda hat) to get a betastar hat, that's a legitimate "trick" to avoid dealing with the parameter being involved in both
coefficients. Bit this is only a valid "trick" when estimating it classically I think.

I don't think there is anything illegal about it. That is what I essentially suggested when I defined the parameter as beta1mrho and made beta a generated quantity. In fact, it is a bit odd in a maximum likelihood framework because the sampling distribution of the ratio would be Cauchy-like but that is not a problem in the full Bayesian setting.

Ben

 

[Mark Leeds]

Thanks Andrew, Bob and Ben for your comments. Going forward, I will refer to bob's formulation as the "wiki" koyck and the formulation with rho on both coefficients as the "depooter" koyck. I don't think
we need to deal with the baltagi version. It only differs from the depooter in the initial assumption about the error term on the long infinite coefficient formulation.

A) Ben:  I will try the "backing out" of betahat method using bob's formulation and see how that goes. Note though that the short form of the model in the wiki is incorrect because there is no transformation that
will take you from the infinite coefficient formulation at the top of the wiki to the reduced form that they have in the wiki. I think the code you sent  is related or the same as the backing out method. I have to look
at it more carefully. I'm just not so convinced it's valid to do that in the bayesian setting because rho being involved in the X_t coefficient, should effect the joint posterior. But I'll worry about understanding
that later.

B) Andrew: Thanks for the ARM centering reference and paper recommendation. I will read the relevant pages in gelman and hill. Also, I'd love to nail this problem and get an applied paper with whomever else was interested. It was a long fight previously that knocked me down for the count and I'd love to get up and end up winning the fight !!!!!!!

C)  Bob:  I have not seen so many centering applications in time series but that doesn't mean it's not used. I'll look around. I will also try centering on the actual depooter model to see if that can help deal with rho being involved in both coefficients. I tend to think it won't though. Also, thanks for your code.

                                                                                                                                                                                                        Mark

                                                                                                                                                                Mark

[Bob Carpenter]

On 4/9/13 12:05 PM, Mark Leeds wrote:
> Thanks Andrew, Bob and Ben for your comments. Going forward, I will refer to bob's formulation as the "wiki" koyck and
> the formulation with rho on both coefficients as the "depooter" koyck. I don't think
> we need to deal with the baltagi version. It only differs from the depooter in the initial assumption about the error
> term on the long infinite coefficient formulation.

Overall, the model won't make sense in a Bayesian setting with the
parameterization (1-lambda)* beta * x[n]. The problem is that
beta needs to get huge when lambda gets close to 1, so what kind
of prior makes sense for it?

Putting a prior on beta makes much more sense.

> A) Ben: I will try the "backing out" of betahat method using bob's formulation and see how that goes. Note though that
> the short form of the model in the wiki is incorrect because there is no transformation that
> will take you from the infinite coefficient formulation at the top of the wiki to the reduced form that they have in the
> wiki. I think the code you sent is related or the same as the backing out method. I have to look
> at it more carefully. I'm just not so convinced it's valid to do that in the bayesian setting because rho being involved
> in the X_t coefficient, should effect the joint posterior. But I'll worry about understanding
> that later.

The issue is that you're sampling here. You can transform samples
and still get valid samples.

For example,

y ~ normal(0,1);
z <- y * sigma + mu;

z will have a normal(mu,sigma) distribution in the samples.

You can't do it the other way around
without accounting for the Jacobian. For instance, this
won't do what you expect in tems of samples for y.

z <- (y - mu) / sigma;
z ~ normal(0,1);

For this, you need to adjust for the Jacobian because the
variable being sampled (on the left of the ~) is transformed.
(It is kosher if sigma is a constant rather than another parameter,
because the log Jacobian determinant is log |1/sigma|, which is
still constant if sigma is a constant.)

In Stan, that adjustment is the following (assuming sigma is
constrained to be positive):

lp__ <- lp__ - log(sigma);

> B) Andrew: Thanks for the ARM centering reference and paper recommendation. I will read the relevant pages in gelman and
> hill. Also, I'd love to nail this problem and get an applied paper with whomever else was interested. It was a long
> fight previously that knocked me down for the count and I'd love to get up and end up winning the fight !!!!!!!
>
> C) Bob: I have not seen so many centering applications in time series but that doesn't mean it's not used. I'll look
> around. I will also try centering on the actual depooter model to see if that can help deal with rho being involved in
> both coefficients. I tend to think it won't though. Also, thanks for your code.

Most of the time series I've seen assume centered data. But
this may be a coincidence due to what people have sent me and
what I've read in the Wikipedia.

For instance, in the paper someone else referred me to by
Kim, Shepherd and Chib on stochastic volatility and GARCH
models:

http://finance.martinsewell.com/stylized-facts/volatility/KimShephardChib1998.pdf

They call it "mean corrected return", but it's the same idea.
It's characterized by the lack of an intercept term in the regression.

Most descriptions of Gaussian process models assume centered data (such
as most of the discussion in the Rasmussen and Williams book).

- Bob

[Mark Leeds]

thanks bob. I'll print out and read carefully but one thing about what you said about it being a valid model ( and you are probably right about centering in time series. I'm not totally familar with
the whole space ).

:the real koyck distributed lag model is actually the long form with the infinite sum. the koyck transformation reduces it to the simpler form with the beta*(1-lambda). So, it has to be valid ( maybe not in the bayesian sense ) econometrically speaking because the koyck distributed lag is actually

A)  y_t = beta * (sum from i = 1 to infinity of lamba^i X_t-i ) + e_t So, a beta coefficient times an infinite sum of past coefficients each with coefficient_i equal to lambda to the i.
 
The reduced form of the koyck model  above is

B) y_t = lambda*y_t-1 + beta*(1-lambda)*X_t.  + v_t ( v_t normal assuming v_t = e_t - lambda e_t-1 )


B) is just the transform of A) using the koyck transformation. You can see the baltagi pdf or adl_gibbs.pdf for how one gets from the long form to the shorter form.  So, it is valid
econometrically. Now the adl_gibbs.pdf paper claims that the model above can be estimated using gibbs sampling. It talks somewhat separately about problems at the boundary and the
gelman-meng paper involving bimodal posteriors and posteriors with ridges etc.

You're understanding of STAN and bayesian is way greater than mine so you may be interested in reading pages 1-23 of adl_gibbs.pdf.

The first 10 or so pages talk about gibbs sampling and about possible problems at the boundary and  also the gelman-meng paper involving bimodal posteriors and posteriors with ridges etc.

Then the next 10 pages or talk about various regression models and how they can be estimated using gibbs sampling. The last model ( 19-22) is the koyck distributed lag I describe above.
They claim that the model can be ( and is ) estimated using gibbs sampling. They get very nice conditional densities and claim to sample from them to estimate
the model. This I doubt because I tried it.

But then, on pages 22-23, they kind of wave their hands and do talk about the identifiability issue when rho is near 1 and how special priors for rho might be able to take care of it etc.  This part is a little
unclear to me and you may be able to understand it a lot better. and it may be possible to implement in STAN ?

So, overall the paper is kind of a little unclear because they claim to estimate the model using gibbs sampling but then in the next breath, they talk about how rho near
1 causes problems. I think the paper was published somewhere ( there's also an error in equation 44. this could have contributed to gibbs sampling working possibly ? ) but like I said,
I doubt the results are correct. On the other hand, if we could understand the "special prior" they recommend using in order to deal with the boundary issue, then maybe it is possible  using STAN?

I have an email somewhere from Lennart Hoogerhiede, author of Admit, who sent me his opinion and is an MCMC guru expert like yourself. 'll dig it up and send you what he said. Hold on a moment. and thanks for your comments. Actually, I'll send it later. I referee basketball at night and I have to leave because I'm late. Thanks again for below.

[Ben Goodrich]

On Tuesday, April 9, 2013 4:45:51 PM UTC-4, Bob Carpenter wrote:

On 4/9/13 12:05 PM, Mark Leeds wrote:
> Thanks Andrew, Bob and Ben for your comments. Going forward, I will refer to bob's formulation as the "wiki" koyck and
> the formulation with rho on both coefficients as the "depooter" koyck. I don't think
> we need to deal with the baltagi version. It only differs from the depooter in the initial assumption about the error
> term on the long infinite coefficient formulation.

Overall, the model won't make sense in a Bayesian setting with the
parameterization (1-lambda)* beta * x[n].  The problem is that
beta needs to get huge when lambda gets close to 1, so what kind
of prior makes sense for it?

Putting a prior on beta makes much more sense.

I think it is actually the opposite: when the autoregressive parameter (lambda) gets close to 1, beta needs to go to zero. Technically, beta isn't identified by the data in that case, but presumably you want the prior to be such that the posterior distribution of beta is near zero so that you conclude the outcome is a random walk.

So, you could specify a normal or t prior on beta where the prior location is zero and the prior scale depends multiplicatively 1 - lambda^2. But I still think it is better to parameterize the coefficient on x as beta1mlambda, which stands for beta * (1 - lambda), again with a normal or t prior on beta1mlambda where the prior location is zero and the prior scale depends multiplicatively 1 - lambda^2. Then make beta a generated quantity.

Ben

[Bob Carpenter]

On 4/9/13 7:32 PM, Ben Goodrich wrote:
> On Tuesday, April 9, 2013 4:45:51 PM UTC-4, Bob Carpenter wrote:
>
> On 4/9/13 12:05 PM, Mark Leeds wrote:
> > Thanks Andrew, Bob and Ben for your comments. Going forward, I will refer to bob's formulation as the "wiki"
> koyck and
> > the formulation with rho on both coefficients as the "depooter" koyck. I don't think
> > we need to deal with the baltagi version. It only differs from the depooter in the initial assumption about the
> error
> > term on the long infinite coefficient formulation.
>
> Overall, the model won't make sense in a Bayesian setting with the
> parameterization (1-lambda)* beta * x[n]. The problem is that
> beta needs to get huge when lambda gets close to 1, so what kind
> of prior makes sense for it?
>
> Putting a prior on beta makes much more sense.
>
>
> I think it is actually the opposite: when the autoregressive parameter (lambda) gets close to 1, beta needs to go to
> zero.

?? The Wikipedia had a simple coefficient in the regression:

beta_raw * x[n]

If we change the model to include the term

beta * (1 - lambda) * x[n]

then

beta = beta_raw / (1 - lambda).

And so if beta_raw > 0, then as lambda -> 1, we have
beta -> +infinity.

Or am I totally missing something here? (Wouldn't be
the first time!)

- Bob

[Ben Goodrich]

If beta_raw were independent of lambda, then yes. But I don't think we should consider these priors to be independent. If you write it like

beta = beta_raw(lambda) / (1 - lambda)

where beta_raw is a function of lambda, and want beta to approach 0 as lambda approaches 1, then making the prior scale be (2.5 * (1 - lambda^2))^2  is one way to do that.

lambda <-seq(from = .Machine$double.eps, to = 1 - .Machine$double.eps, length.out = 1000)
foo <- function(l) sd(rnorm(10000, 0, (2.5 * (1 - l^2))^2) / (1 - l))
plot(lambda, sapply(lambda, foo), type = "l", ylab = "SD of beta") # attached

Ben

[Carlos]

Ok, sorry for getting back sooo late to this. This might be stupid, but I think that the problem might be the way that Mark is generating the artifitial data. He gives the X series as a sequence of integers that goes from 1 to N. This is just imposing a trend on Y, and therefore, the model will not be able to be estimated probably unless he stationarizes the time series, which is the point of the model. I attach here other version of the code where X is generated just as an iid, and even for a close to unit root rho, it works reasonably well...

[Carlos]

PS: Yes, in time series, people usually work with zero-mean data (we stationarize the data before doing anything else)

[Mark Leeds]

Hi Bob: I don't want to comment on beta_raw because beta_raw is from a model that can't be deduced from the original koyck model. ( the long form where there are an infinite number of coefficients ).
The model in the wiki is not possible in the sense of being a model that can be obtained from a transformation of the original infinite sum formulation. That's important because the whole point of obtaining
the reduced form of koyck distributed lag is to be able to estimate the parameters in THE ORIGINAL LONG FORM OF THE MODEL. In other words, the long form is a model that captures dampening behavior where the weights of the coefficients get less and less weight the further you go back in time. But estimating that model is difficult directly. So, the koyck transformation results in the "short form" that is
much easier to estimate ( atleast classically. but you can run into problems there too if you don't use numerical optimization and gradients ).  The point I'm trying to make is that the short form of model in the wiki doesn't have a long form correspondence. It's incorrect in terms of it corresponding to anything useful from a "distributed lag" viewpoint.  Wiki's can definitely be wrong. I've seen this happen before.

So, getting back to "Depooter" koyck:

The simplistic way I see it is that   lambda going to one in the "depooter" koyck causes an "indentifiability" problem because beta*(1-lambda) gets close to zero so it's impossible to estimate beta.

but let me write here what Lennart wrote to me in an email.  I don't have a hardcopy of the email so I'll just write it. and Lennart definitely knows what he's doing so I would bet there's some validity to the following:

Quoting Lennart from May 16, 2010 in an email to Mark Leeds:

" Hi Mark: I read your email and the relevant pages in the paper by Depooter et al: I think it is the prior that causes the problems. Indeed there must be a ridge in the likelihood for rho equal to 1 which I think is
not integrable ( having integral of infinity ),  if the prior does not "correct" it. So, under the uniform prior for rho,  the posterior kernel does not correspond to a proper distribution. In other words, intuitively
stated, the "scaling  constant" should be infinity which ruins the results.

I think a regularization prior for rho such as probability(rho) ~ (1-rho hat)^2 which keeps rho somewhat away from 1, ( similar to equation (52) in Depooter et al ) would solve the problems.".

End of Quote from Lennart Hoogerheide.

I can't remember what happened after this in terms of whether I tried this with the Admit package or not. I don't think I did because I didn't understand how that could be a prior because I thought
the prior needed to integrate to 1. At that point, I think I decided to give up and go the classical route using numerical optimization and gradients.

Anyway, maybe above sheds some light on the issue. It's pretty similar to what's talked about on pages 22-23 of adl_gibbs.pdf.

Ben: I'll read your email tomorrow. It's late and my brain hurts.

[Bob Carpenter]

On 4/10/13 1:54 AM, Mark Leeds wrote:
>
> Hi Bob: I don't want to comment on beta_raw because beta_raw is from a model that can't be deduced from the original
> koyck model. ( the long form where there are an infinite number of coefficients ).

It's just a reparameterization of what you had, which
will help with arithmetic stability in the estimate.

Given what you say below, it shouldn't be necessary.

> The model in the wiki is not possible in the sense of being a model that can be obtained from a transformation of the
> original infinite sum formulation. That's important because the whole point of obtaining
> the reduced form of koyck distributed lag is to be able to estimate the parameters in THE ORIGINAL LONG FORM OF THE
> MODEL. In other words, the long form is a model that captures dampening behavior where the weights of the coefficients
> get less and less weight the further you go back in time. But estimating that model is difficult directly. So, the koyck
> transformation results in the "short form" that is
> much easier to estimate ( atleast classically. but you can run into problems there too if you don't use numerical
> optimization and gradients ).

For maximum likelihood, there shouldn't be an issue until
lambda gets so close to 1 that beta overflows to keep up.

> The point I'm trying to make is that the short form of model in the wiki doesn't have a
> long form correspondence. It's incorrect in terms of it corresponding to anything useful from a "distributed lag"
> viewpoint.

I'm not even sure what people use these models for!
I'm just commenting on the Stan formulation and the
behavior you can expect from different parameterizations
of your model.

> Wiki's can definitely be wrong. I've seen this happen before.

Of course --- that's why I was suggesting changing
the Wikipedia entry!

> So, getting back to "Depooter" koyck:
>
> The simplistic way I see it is that lambda going to one in the "depooter" koyck causes an "indentifiability" problem
> because beta*(1-lambda) gets close to zero so it's impossible to estimate beta.

Why does this make it impossible to estimate beta?
If there's an effect for the predictor, then beta should
just grow to compensate.

A Bayesian prior on beta (such as beta ~ normal(0,1000))
will stop beta from growing enough, so the
beta * (1 - rho) term will just go to zero as
lambda goes to 1.

> but let me write here what Lennart wrote to me in an email. I don't have a hardcopy of the email so I'll just write it.
> and Lennart definitely knows what he's doing so I would bet there's some validity to the following:
>
> Quoting Lennart from May 16, 2010 in an email to Mark Leeds:
>
> " Hi Mark: I read your email and the relevant pages in the paper by Depooter et al: I think it is the prior that causes
> the problems. Indeed there must be a ridge in the likelihood for rho equal to 1 which I think is
> not integrable ( having integral of infinity ), if the prior does not "correct" it. So, under the uniform prior for
> rho, the posterior kernel does not correspond to a proper distribution. In other words, intuitively
> stated, the "scaling constant" should be infinity which ruins the results.

> I think a regularization prior for rho such as probability(rho) ~ (1-rho hat)^2 which keeps rho somewhat away from 1, (
> similar to equation (52) in Depooter et al ) would solve the problems.".

:-) I thought you wanted rho to go to 1 if the data
supported it. If you really want to impose a strong enough
prior to bound rho away from 1, that will manage the scale
of beta (to the degree that rho is bound away from 1).

> End of Quote from Lennart Hoogerheide.
>
> I can't remember what happened after this in terms of whether I tried this with the Admit package or not. I don't think
> I did because I didn't understand how that could be a prior because I thought
> the prior needed to integrate to 1.

The integral just needs to be finite, which is equivalent to
integrating to 1, but dropping the proportionality constant.
In this case it's easy because rho is constrained to [0,1].

> At that point, I think I decided to give up and go the classical route using
> numerical optimization and gradients.
>
> Anyway, maybe above sheds some light on the issue. It's pretty similar to what's talked about on pages 22-23 of
> adl_gibbs.pdf.
>
> Ben: I'll read your email tomorrow. It's late and my brain hurts.

I think after you read Ben's mail you'll see we're all on
the same page again. Ben also had some suggestions on the
non-independence of the prior on rho and beta. You could
also overcome the issue that way, where the scale of beta is
allowed to grow as rho approaches 1.

- Bob

[Mark Leeds]

Thanks Bob. I'm doing a sports marathon today ( softball umpiring then basketball  refereeing) so I can't read below carefully now. But I did see that we're all on the same page part. So that's good. So, let
me try whatever you're suggesting ( I think it's the "Ben" betamrho approach. !!!!!! )  and I'll let you know how it goes. I still haven't even tried it ( nor read Ben's emails ). Maybe that will be fine and will close the deal.  I'm really anxious to look at this today but today there are some sports committments that I gotta attend to. Hopefully late tonight.

Also, I didn't understand what you meant by using betamrho ( I didn't know doing that was valid in the sense of backing out betastar ) so that's why I was trying to stress that the wikl was wrong.
I thought you were actually trying to estimate the wiki model. I'll also explain what the  depooter koyck model is used for another time. In mostly for advertising/sales stuff.

Finally, to answer your question: I want rho to be estimated correctly at whatever the true value is. ( near one, not near one, whatever ). But what happens with my STAN formulation currently is that, as the true value of rho goes near 1.0 ( actually above 0.75 ), the beta estimate gets whackey. Thanks to you and Ben for your help and have a good day and I'll get back to both of you hopefully tomorrow. Carlos,
I didn't respond to you either but I'm pretty certain that there's nothing wrong with my simulation. For the wiki model, Ben's code worked perfecto and the artificial simulation is the same as mine for horseshoe
and hand grenade purposes. But I'll read your email also and get back. I only saw the very beginning.

[Mark Leeds]

I have to leave again so I'll make this brief while hopefully still being informative. the R code and STAN code associated with this recent experiment  is at the bottom of the email.

A) I really just implemented and ran what Ben suggested originally but, if someone could check my code, it's appreciated because I haven't read the manual fully so I don't know yet when order matters and when it doesn't. Amazing manual by the way.  Whoever wrote that should e-mail R-core and ask them if they'll pay he/she to document R similarly !!!!

B) I ran a similar scenario set to what I ran earlier in terms of sets of triplets of parameters ( each triplet only varies in beta values )  and then as I go down, I keep decreasing lambda to see how that effects
the beta estimates.

C) The results are better. The estimates of beta don't get totally whacky like they did in the earlier experiment but for lambda ( which of course used to be rho. we just changed notation somewhere along the way ) above 0.9, the beta estimates are still not very good. I put the beta estimate at the top of each triplet and gave it a "review" so that it's not necessary to re-run all the code. I set.seed(123) at the top
if the code before running each scenario but there's some randomness of course going on in STAN so the results won't be the same. It seems that below lambda = 0.9, the estimates pretty much
remain as good as they're going to get for the most part. In my other run where I explicitly code lambda(1-rho), the improvement seemed to be less spikey and more gradual. Not sure if that means
anything.

D) I don't claim to be an HMC expert but I'm trying to think about something that my understanding is not 100 percent clear on so hear me out if you don't mind. Lennart says that the log posterior goes to
infinity as rho goes to 1. I'm not clear on how he knows that but, assuming it's true: could the problem then correspond to the example on pages 30-31 of "Probabilistic Inference Using Markov Chain Monte Carlo
Methods" by Radford Neal. If someone could take a look at that discussion and the picture in c). The density is spiking up in one place and the values of A in that place are all very low so MCMC sampling will
be problematic. I don''t know enough about HMC ( atleast yet. as time permits, I'm going back over then DUANE paper for the 38th time ) yet to say anything with confidence but maybe the problem has something to do with that which I think is what Lennart was saying ?

So, given D) is true ( which it may not be ), I have two thoughts/ideas:

1) restrict rho to be uniform(0,0.95). only problem with that is you leave out cases where rho really is greater than that.

2) use the prior that lennart suggested to force the posterior density to be zero when rho gets near 1.0.

Thanks for any other comments/suggestions/wisdom. I'm still getting my brain immersed in HMC  yes, the trajectory stopping thing is currently beyond me but I still have to go back and see where the actual
posterior density is involved in HMC. I think it's only involved because it's the gradient of the posterior that's in the differential equations that simulate states of parameters. but if that's the case,  that's so different from regular metropolis where the actual posterior atleast gets plugged into on the top and bottom in the acceptance step. The metropolis argument makes perfect sense to me.  In HMC, the metrpolis type
acceptance is used  the actual posterior isn't even involved in the acceptance step !!!!!! and this continues to confuse me in terms of how one can just use multivariate normals ( I understand that the normals are conditional but that the conditionals are just the marginals because of independence of the two energies in the expression for the hamiltonian. I just don't understand why one don't plug into the posterior
like would be done in a regular metropolis step. that has something to do with sampling from the differential equations rather than from the density itself but I'm working on it  )  but I'm getting it slowly. 

So my point is that D)  is all guess work.I'll be out the rest of the day so no immediate replies. And no should try to explain HMC to me. I keep reading and reading and reading and when  I get to the point where I can't understand it any better, I'll ask someone on the list. Not that I wouldn't appreciate an explanation but I'd rather get one when I'm done making my own attempts.  Thanks for reading this and any comments about A-D and 1) and 2).

# R CODE FOR BEN'S SUGGESTION
#==================================================================================================

library("rstan");

set.seed(123)

T <- 200;
x <- rnorm(T,0,1);

#=========================================

# NOT OKAY : 0.2

#=========================================
#beta <- 1.0;
#lambda <- 0.95
#sigma <- 1.0

#=========================================
# NOT OKAY : 1.5
#=========================================

#beta <- 2.6;
#lambda <- 0.95
#sigma <- 1.0

#=========================================
# NOT OKAY : -0.1
#=========================================

#beta <- 0.5;
#lambda <- 0.95
#sigma <- 1.0
#====================================

# NOT OKAY: 0.6

#====================================
#beta <- 1.0;
#lambda <- 0.90
#sigma <- 1.0

#====================================
# NOT GREAT : 1.8
#========================================

#beta <- 2.6;
#lambda <- 0.90
#sigma <- 1.0

#========================================
# NOT SO GOOD : 0.2
#========================================

#beta <- 0.5;
#lambda <- 0.90
#sigma <- 1.0
#========================================

#NOT BAD: 0.8

#=========================================
#beta <- 1.0;
#lambda <- 0.85
#sigma <- 1.0

#=========================================
#SO SO : 2.0
#=========================================

#beta <- 2.6;
#lambda <- 0.85
#sigma <- 1.0

#=========================================
#SO SO : 0.3
#=========================================

#beta <- 0.5;
#lambda <- 0.85
#sigma <- 1.0
#====================================

# SO SO:  0.7

#====================================
#beta <- 1.0;
#lambda <- 0.80
#sigma <- 1.0

#====================================
# DECENT: 2.2
#====================================

#beta <- 2.6;
#lambda <- 0.80
#sigma <- 1.0

#====================================
#SO SO : 0.3
#====================================

#beta <- 0.5;
#lambda <- 0.80
#sigma <- 1.0
#====================================

# OKAY:  0.8

#====================================
#beta <- 1.0;
#lambda <- 0.75
#sigma <- 1.0

#====================================
# GOOD: 2.7
#====================================

#beta <- 2.6;
#lambda <- 0.75
#sigma <- 1.

#====================================
# STILL SO SO : 0.3
#====================================

#beta <- 0.5;
#lambda <- 0.75
#sigma <- 1.0
#====================================

# OKAY:  0.8
#====================================
#beta <- 1.0;
#lambda <- 0.70
#sigma <- 1.0
#====================================
# GOOD: 2.4
#====================================

#beta <- 2.6;
#lambda <- 0.70
#sigma <- 1.0

#====================================
# GOOD: 0.4
#====================================

#beta <- 0.5;
#lambda <- 0.70
#sigma <- 1.0

#===========================================


y <- rep(NA,T);

y[1] <- rnorm(1,beta * (1 - lambda) * x[1],sigma);

for (t in 2:T)
  y[t] <- rnorm(1, beta * (1 - lambda) * x[t] + lambda * y[t-1], sigma);

fit <- stan(file="koyck_gen.stan", data=c("T","y","x"), iter=1000, chains=4);
print(fit)

#=========================================================================================
# STAN CODE
#===============================================================================================================================================

data {
  int<lower=0> T; // number of time points
  real y[T]; // output at time t
  real x[T]; // predictor for time t
}
parameters {

  real beta1mlambda; // equals beta * (1-lambda)

  real <lower=0, upper=1> lambda; // lag
  real <lower=0> sigma; // noise scale
}

model {

  beta1mlambda ~ cauchy(0,5);

  lambda ~ uniform(0,1);
  sigma ~ cauchy(0,5);

  y[1] ~  normal(beta1mlambda * x[1], sigma);

  for (t in 2:T)

    y[t] ~ normal( beta1mlambda * x[t] + lambda * y[t-1],sigma);
}

generated quantities {
  real beta;
  beta <-  beta1mlambda / (1 - lambda);

[Michael Betancourt]

> D) I don't claim to be an HMC expert but I'm trying to think about something that my understanding is not 100 percent clear on so hear me out if you don't mind. Lennart says that the log posterior goes to
> infinity as rho goes to 1. I'm not clear on how he knows that but, assuming it's true: could the problem then correspond to the example on pages 30-31 of "Probabilistic Inference Using Markov Chain Monte Carlo
> Methods" by Radford Neal. If someone could take a look at that discussion and the picture in c). The density is spiking up in one place and the values of A in that place are all very low so MCMC sampling will
> be problematic. I don''t know enough about HMC ( atleast yet. as time permits, I'm going back over then DUANE paper for the 38th time ) yet to say anything with confidence but maybe the problem has something to do with that which I think is what Lennart was saying ?

It sounds like you're dramatically over-thinking HMC, so let me try to explain.

In HMC you're not sampling from your target distribution pi(q), but rather an extended distribution pi(p, q) = pi(p | q) p (q). Hamiltonian dynamics provides a Metropolis proposal for {p, q} which is accepted with the usual Metropolis acceptance procedure for the extended distribution pi(p, q).

You then get a sample from pi(q) by simply ignoring the variables p since the marginal of the extended distribution is the target distribution, \int dp pi(p, q) = pi(q) \int dp pi(p | q) = pi(q).

I haven't looked at your model in any depth, but if there is finite posterior density at rho = 1 then there will be problems sampling with HMC. The problem is that when we transform rho to an unconstrained space the boundaries go to infinity, [0, 1] -> [-inf, inf]. If the posterior density is finite at rho = 1 then it will also be finite at rho_unconstrained = inf, and Hamiltonian dynamics falls apart when the Hamiltonian is non-zero at infinity. Computationally it's easy to see why -- if the density is non-zero out there, then the chain needs to travel all the way out to infinity which would take an infinite amount of time, not to mention wreaking havoc on any numerical implementation. This is case where using bouncing to constrain variables can help, but we haven't implemented that in Stan.

-Mike

[Bob Carpenter]

On 4/11/13 3:53 PM, Michael Betancourt wrote:
>> D) I don't claim to be an HMC expert but I'm trying to think about something that my understanding is not 100 percent clear on so hear me out if you don't mind. Lennart says that the log posterior goes to
>> infinity as rho goes to 1. I'm not clear on how he knows that

I'll guess calculus. Just calculate

lim_{rho -> 1} log p(rho,...)

> but, assuming it's true:

If you can't verify this yourself, you'll have a hard
time understanding HMC.

- Bob

[Bob Carpenter]

On 4/11/13 2:59 PM, Mark Leeds wrote:
> I have to leave again so I'll make this brief while hopefully still being informative. the R code and STAN code
> associated with this recent experiment is at the bottom of the email.
>
> A) I really just implemented and ran what Ben suggested originally but, if someone could check my code, it's appreciated
> because I haven't read the manual fully so I don't know yet when order matters and when it doesn't. Amazing manual by
> the way. Whoever wrote that should e-mail R-core and ask them if they'll pay he/she to document R similarly !!!!

We're happy to answer questions, but please read
the manual first. We don't have time to recapitulate
the manual on the mailing list for every user.

- Bob

[Ben Goodrich]

On Thursday, April 11, 2013 2:59:13 PM UTC-4, Mark Leeds wrote:

#=========================================================================================
# STAN CODE
#===============================================================================================================================================


data {
  int<lower=0> T; // number of time points
  real y[T]; // output at time t
  real x[T]; // predictor for time t
}
parameters {
  real beta1mlambda; // equals beta * (1-lambda)
  real <lower=0, upper=1> lambda; // lag
  real <lower=0> sigma; // noise scale
}

model {
  beta1mlambda ~ cauchy(0,5);
  lambda ~ uniform(0,1);
  sigma ~ cauchy(0,5);
  y[1] ~  normal(beta1mlambda * x[1], sigma);
  for (t in 2:T)
    y[t] ~ normal( beta1mlambda * x[t] + lambda * y[t-1],sigma);
}

generated quantities {
  real beta;
  beta <-  beta1mlambda / (1 - lambda);
}

I think in the PDF you sent earlier, the authors are presuming an improper prior on beta. In that case, as lambda -> 1, the marginal posterior distribution of beta is not well-defined. However, with a proper prior on beta (or beta1mlambda), then the marginal posterior distribution is just the prior distribution in the case where lambda -> 1. That would be true for Bayesian inference in general; it doesn't have anything in particular to do with HMC or Gibbs sampling or whatever algorithm you use to draw from the posterior.

I think most people would want the posterior distribution of beta to converge to 0 as lambda -> 1, which implies that the prior for beta1mlambda can't be independent of lambda. I would do this:

data {
  int<lower=0> T; // number of time points
  real y[T]; // output at time t
  real x[T]; // predictor for time t
}
parameters {
  real beta1mlambda; // equals beta * (1-lambda)
  real <lower=0, upper=1> lambda; // lag
  real <lower=0> sigma; // noise scale
}

model {

  real mu[T];
  beta1mlambda ~ normal(0,square(2.236 * (1 - square(lambda))));
  // implicit: lambda ~ uniform(0,1);
  sigma ~ cauchy(0,5);
  mu[1] <- y[1]

  for (t in 2:T)

    mu[t] <- beta1mlambda * x[t] + lambda * y[t-1];

  y ~ normal(mu, sigma);

}
generated quantities {
  real beta;
  beta <-  beta1mlambda / (1 - lambda);
}

Ben

[Mark Leeds]

thanks michael. printing out and will read carefully. thanks.

[Mark Leeds]

thanks bob. I know how to take limits and logs but I don't know what the joint posterior density is , so I need that function to take  the log and the limit. I'll look at the depooter paper again because there
is some material in there on 22-23 about the hessian of the likelihood. maybe that will clear up my confusion. but I was never even clear how they got their joint prior. So the joint posterior probably won't be that
much clearer. I'll take a look.

[Mark Leeds]

okay bob. my apologies for consuming your time. I'll make sure I read the manual front to back before I ask another stan language question.

[Mark Leeds]

[Mark Leeds]

thanks ben. I'll take a look at the adl_gibbs.pdf again but I remember being confused by the joint prior. I've read the paper many times but not at all recently. Thanks for your suggestions and code below. I'll get back to you. Also, do you have an email that I can communicate directly ? I promise not to over-burden you and it would be better to not burden others on the list with my confusion/stupid questions. thanks.

[Bob Carpenter]

The posterior density is proportional to the joint
density by Bayes's rule. So it's just

lim_{rho -> 1} p(rho,...) * p(data|rho,...)

where ... are the other parameters.

- Bob

[Mark Leeds]

thanks bob but I do know bayes rule. I've been telling everyone to read adl_gibbs.pdf when I've been guilty myself of not having read it in 3 years. So I'll check it out to see what p(data|rho) is and report
back with any useful findings. I do know that it blows up as lim rho -> 1.

[Mark Leeds]

Hi Carlos: No, my simulation code doesn't have to do with adding a trend. That was just my silly deterministic way of generating random coefficients to test out how well the parameters are estimated.
It can be confusing because of the t-subscripting but the reduced form of the koyck model ( the one being coded in STAN )  is really just a  cross sectional regression. ( although there are issues with
correlation between regressor and noise term in estimating it as a regression but I shouldn't go into that here ). I write the model for each time t below:

y_1 = lambda* y_0 + lambda((1-p)* X_1 + epsilon_1
y_2 = lambda*y_1 + lambda(1-p)* X_2 + epsilon_2
y_3 = lambda * y_2 + lambda*(1-p)*X_3  + epsilon_3
.
.
.

the subscripts are t ( as in time ) but you can really just think of them as indices. So, I just incremented the X_t's by 1 each time in the code but that's not a trend model as you see from above. Just
different X coefficients in order to see how the parameters get estimated. If there was a second index in the y vector that was indexed by t, then that would be assuming some kind of trend. 
If you want to think of above as a time series model,  you actually can think of it as an AR(1) with an exogenous coefficient but that's not really a useful way to think of it atleast in my opinion.
You can shuffle around the (X,y ) pairs in the code using "sample" in R and re-run the model and you'll see that you'll get the same estimates. I'm not sure if that explains your trend question. Let me know.

But, the true model ( what I refer to as the long form that is really the true model. short form is just a reparameterization in order to make estimation easier  ) underlying the short form above is actually a
lagged time series model but still with no time trend. Specifically:

y_t = beta*(sum from i = zero to infinity of (lambda^i *X_t-i) + epsilon_t  for every t

This is the really interesting formulation of the model and the way I like to think of it. ( I also sometimes think of it in the short form way but this one makes it more clear to me what's really happening ).

Note that It's the same model as the short form but re-written. This model captures ( or is supposed to capture anyway ) the lagged effects of X_t  on y_t over time. So, for example, say you had a shock X_t at time t where X_t was an abundance of dollars spent on advertising at time t. ( company is doing a special advertising push say at time t ). Also, suppose you let y_t denotes sales at time t. Then what this model tries to do is capture the lagged effect of this $ investment as time goes forward from time t. Since the weights are decaying as you go firhter back in time, it's assumes that as you go back further in time,
the effect on y_t decreases. So, it's kind of like a decreasing momentum model. But it's neat because most econometric/time series models don't handle effects that occur over time. They assume the effect is all at once.

There are a lot of details that I'm leaving out but there are tons of papers on the net about this model. If  you don't know which one to read, let me know and I can dig up something that's good and send you the title. I haven't looked at koyck papers about it in a good while but I've read many. The Koyck model is an econometric model that is I think top 5 or top 10 all time in econometric citations.  I hope this helped clear up what's going on in the code and the model.

                                                                                                                                                                                                  Mark













The true use of the reduced form of the koyck model is to construct an estimable model with 3 parameters because, in its long form ( the one of the form beta * sum of infinite lags with each lag multiplied by weight lambda^i ), is extremely difficult to estimate. So, koyck figured the transformation out which leads to the reduced version which is what we're trying to estimate because it's also 3 parameters
but much easier to estimate.

 

On Tuesday, April 9, 2013 11:19:05 PM UTC-4, Carlos wrote:

[Mark Leeds]

correction on typo below. It should read:

y_1 = lambda* y_0 + beta((1-lambda)* X_1 + epsilon_1
y_2 = lambda*y_1 + beta(1-lambda)* X_2 + epsilon_2
y_3 = lambda * y_2 + beta*(1-lambda)*X_3  + epsilon_3

I confused my lambdas, rhos and betas !!!! my apologies. i will trash rho and use lambda going forward.

[Mark Leeds]

Hi Michael : Below was a beautiful explanation. I think I'm getting closer and you're also correct that I was over thinking it or atleast wrongly thinking it !!!!. If you can bear with me and
just explain this last piece at your convenience, I'd appreciate it. ( this is not the trajectory stuff but that's a future battle ).  So, here goes. I will state what I think I understand first and then where my confusion
is. But, if I'm incorrect in where I think I'm correct, then please let me know that also of course.

In regular metropolis hastings, ( local, global, whatever ),  the proposal for the new state is generated using some proposal distribution. ( I'm aware of the difficulties of which one to choose and the variance issue etc ). Then, one  uses the acceptance step as a way to insure detailed balance. Finally, detailed balance implies convergence of the constructed markov chain.to the stationary distribution. Because of the way the acceptance prob is constructed, min(1, prob(new state)prob(old state|new state) divided by prob(old state)prob(new state|old state) ), this stationary distribution fortunately corresponds to the one we are trying to generate from in the first  place !!!! Good,  got that.

Now, in HMC, ( the last time I studied phyiscs was 1984 and there were no differential equations so there's going to be hand waving going on here ) we are trying to draw from the distribution of q which I denote
as p(q) but we actually have to two sets of variables, q and p, where the p are the artificial ones just introduced to make HMC work. We can stick to 1 dimensional p and q for what I'm asking here.

Now,  in HMC, the new proposal ( the next p and q)  is generated using the leapfrog steps which are a discretization of the dynamics implied by hamilton's equations. ( except the momentum is actually
generated randomly afterwards to get out of the contour of constant probability but I'll talk about probabilities later ) .  Okay, fine. Then, aside from the way the proposal is done, the acceptance step is exactly the same way it is in metropolis hastings : prob(acceptance of new state) = min(1, prob(new state)prob(old state|new state) divided by [ prob(old state)prob(new state|old state) ].

Then, using the same detailed balance balance argument, the acceptance step implies convergence of the constructed markov chain to the stationary distribution prob(p,q) = prob(p)prob(q|p). But independence
of p and q implies that you can just generate the joint sample of (p, q) and use the generated sample from q because that will converge to the marginal of q. Good. Got that.

But here's where I get lost. According to what I've read,  Prob(p,q) = e to the -[U(q) +K(p)]. ( picture lots of waving of hands here !!!! ).  Now, in the actual statistical problem, we are ultimately trying to sample from some (posterior) density p(q). Like in  the koyck model, we are trying to sample from p(rho, beta, sigma| data) so q = (rho,beta, sigma ).

 So in the actual HMC implementation,  one sets U(q) = - log(prob(rho, beta, sigma).  In this way since prob(p,q) = e to the -(U(q) +K(p)),  this reduces to prob(p,q) = p(q) * [e to the -K(p).]. So ( because of the independence of p and q,  as stated earlier),  the procedure ends up generating samples from the marginal,  prob(q) = prob(rho, beta, sigma) which is exactly what we wanted !!!! Good. Sort of got that,

But, when I read Radford Neal's material, I must be missing some kind of key relationship between the hamilton equations/leapfrog steps and the probabilities. How can one just out
of the seemingly  blue say that Prob(p,q) = e to the -[U(q) +K(p)]. Is that some kind of physical relationship implied by the hamiltonian equations themselves ? In other words, my question is,
do the hamiltonian equations somehow make that relation be true ? that prob(p,q) = e to the -[U(q) + K(p)]. If so, then I'll just accept it ( no pun intended ). I ask this because I don't see it all.
In statistics, one will just out of the blue say that such and such variable X is exponentially distributed or gamma distributed so that one can say that prob(X = x) = some density function. But then this has
to be tested for it's correctness. In HMC, one is stating a probability relationship of variables p and q that are being generated by a set of differential equations. So, I don't think it's an assumption
but rather some kind of by product of the equations ? Is that correct that it's just a fact ? This is my main confusion. if I can understand this, then I think I get it.

Then, to close the story, if it is factual, then is it the discretization of the hamiltonian equations that makes one need the acceptance step. I read somewhere that, if you didn't have to discretize, then you
wouldn't even need the acceptance step because detailed balance would hold exactly without the acceptance step. Again that has to do with the probabilities so it's probably related to the question above.

I guess to summarize I can say that I don't see how one can discretize a differential equation and then have that probability relation be true ? In fact, forget about discretizing. I don't see it even without discretization. But that's where my NULL physics background may be coming into play. Lastly, it's really nice of you to go out of your way to help me. Hopefully you'll be helping some others also ? Thank
you very much.


                                                                                                                                                                                                                            Mark






















Is this somehow implied by the hamiltonian equations ?

On Thursday, April 11, 2013 3:53:24 PM UTC-4, Michael Betancourt wrote:

[Bob Carpenter]

On 4/13/13 8:18 AM, Mark Leeds wrote:
...

> But, when I read Radford Neal's material, I must be missing some kind of key relationship between the hamilton
> equations/leapfrog steps and the probabilities. How can one just out
> of the seemingly blue say that Prob(p,q) = e to the -[U(q) +K(p)]. Is that some kind of physical relationship implied
> by the hamiltonian equations themselves ?

p(p,q) is a density that was defined (nor derived) by
the inventors of HMC. Among other properties, the
definition has these pleasing features:

1. The marginal density p(q) of the position q is the
distribution we want to sample from (by construction).
So if we can sample (p,q) from p(p,q), we can sample q
from p(q).

2. The joint density p(p,q) decomposes into a separable
Hamiltonian form.

In other words, my question is,
> do the hamiltonian equations somehow make that relation be true ? that prob(p,q) = e to the -[U(q) + K(p)]. If so, then
> I'll just accept it ( no pun intended ). I ask this because I don't see it all.

> In statistics, one will just out of the blue say that such and such variable X is exponentially distributed or gamma
> distributed

These are almost always modeling assumptions, which need
to be checked.

> ... is it the discretization of the hamiltonian equations that makes one
> need the acceptance step.

Yes. If we could do the integral corresponding to the Hamiltonian
dynamics perfectly, the acceptance rate would be 100%.

The key to understanding why it maintains detailed balance is due
to volume preservation (this is the part that was tricky for me
to understand).

- Bob

[Mark Leeds]

thanks bob. I'll print out below and read carefully. it's much appreciated. eventually I'll get back to playing with stan/koyck and report back ( will try Ben's idea also ) but I want to understand HMC
better first. have a good day.

[Mark Leeds]

bob: below was extremely helpful particularly when you said "DEFINED" by the inventors of HMC. I get it now I think: The construction of the hamiltionian form for the probabilities is really just a useful way to
say that the probabilities of the two variables ( p and q ) are independent. And, since the hamiltonian form ends up reducing to the density of the one we're interested in TIMES the density of the one we're
not interested in , it's irrelevant if the density of the second one is not true because we don't care about it anyway  !!!!!!!! really neat stuff. The leapfrog steps  are a construction for proposing from the joint of p and q but are connected to the hamiltonian construction through the gradients of the respective densities and therefore "consistent" with the hamiltonian.  WOW.  It feels good when the light goes on !!!!!

Thanks so much to everyone on this list for their patience and help, particularly , Michael, Bob and Ben. I will get back with STAN results but it may be a while because I have some other things currently going on.

[Mark Leeds]

hi all: just for future readers of this thread who are trying to understand HMC like myself, in below I should have said: "the construction of the hamiltonian form for the probabilities is a neat way of enforcing the independence of the density for p and the density for q.

[Mark Leeds]

To all who are following along in the quest for koyck model estimation to allow for rho to be near 1, attached is an earlier version of adl_gibbs.pdf that provides al ofl the gory details as far as  what happens to
the posterior of (lambda,beta) density as lambda approaches 1. If one sits down with pencil and paper, it actually ends up making sense because the posterior derivation involves only integration and a lot of algebra. The paper also explains a fix for the problem so I plan to use exactly their priors along with the fix in STAN and see what happens as rho goes near 1.0. I will also try Ben's earlier recommendation. If either of the latter work, then I will go back to JAGS and see if the same fix works there also. It would be more interesting of course if STAN succeeded and JAGS didn't but we shall see.

                                                                                                                                                       
                                                                                                                                                                                                                    Mark