State space form of time varying AR(1) I would like to implement the model proposed in Dynamic modeling of mean-reverting spreads (Kostas Triantafyllopoulos, Giovanni Montana).
They propose to model a time serie Y_t with the following equations:
(1) Y_t = A_t + B_t * Y_(t-1) + e_t  
(2) A_t = Phi1 * A_(t-1) + nu1_t  
(3) B_t = Phi2 * B_(t-1) + nu2_t

That can be expressed in a state space form:  
(1') Y_t = F_t * theta_t + e_t  
(2') theta_t = Phi * theta_(t-1) + nu_t

with   
F_t = (1, Y_t-1)  
Phi = diag(Phi1, Phi2)

I would like to use R to perform a bayesian update of this model. I have studied the package 'dlm' but in the book Dynamic linear models with R, written by Giovanni Petris (author of the 'dlm' package), it is written (page 113) 

"The matrix F_t of a DLM cannot depend on past values of the observations".   

However, it seams that it is the case in the model proposed above. 
Can someone understand this last sentence and eventually help me perform this implementation in R?
Thank you
Fred
 A: This is a very unusual state space model because the dynamics are included in both the observation equation (1') and the state equation (2'). Usually, the dynamics are only in the state equation and the observation equation is a linear function of the state vector. I don't think any of the state space implementations in R will allow dynamics in the observation equation.
It is possible to re-write the model so that the dynamics are all in the state equation, but then it becomes non-linear.
I suspect you will have to write your own code, or ask the authors if they can give you theirs.
A: You can build this model with AD Model Builder's random effects package.
This is free software available at http://admb-project.org.  What you will
get is full information maximum likelihood solutions with the ability to try
MCMC methods afterwards if you wish. The idea is to regard this as a random
effects problem and integrate over the random effects via the Laplace approximation.  The trick is parameterize if properly so that the Hessian
with respect to the random effects is sparse.  I built the model and a simulator. It seems to work well.  To give you an idea I have included the
ADMB source for the model.
DATA_SECTION
  init_int nobs
  init_vector Y(1,nobs)
  vector resids(2,nobs)
PARAMETER_SECTION
  init_number Phi1
  init_number Phi2
  init_bounded_number log_sigma(-5.0,5.0,2);
  init_bounded_number log_p1(-5.0,5.0,2);
  init_bounded_number log_p2(-5.0,5.0,2);
  objective_function_value f
  init_number A2
  init_number B2
  random_effects_vector A(3,nobs)
  random_effects_vector B(3,nobs)
  vector nu_1(3,nobs);
  vector nu_2(3,nobs);
  vector pred_Y(2,nobs);
  sdreport_number sigma
  sdreport_number p1
  sdreport_number p2
PROCEDURE_SECTION

  f0(log_sigma,log_p1,log_p2,A2,B2);

  f2(3,log_sigma,log_p1,log_p2,A(3),B(3),A2,B2,Phi1,Phi2);

  for (int i=4;i<=nobs;i++)
  {
    f2(i,log_sigma,log_p1,log_p2,A(i),B(i),A(i-1),B(i-1),Phi1,Phi2);
  }

  if (sd_phase())
  {
    sigma=exp(log_sigma);
    p1=exp(log_p1);
    p2=exp(log_p2);
  }


SEPARABLE_FUNCTION void f0( const prevariable& log_sigma, const prevariable& log_p1, const prevariable& log_p2, const prevariable& A2, const prevariable& B2)

  dvariable sigma=exp(log_sigma);
  dvariable p1=exp(log_p1);
  dvariable p2=exp(log_p2);
  f+=square(A2)+square(B2);
  resids(2)=value(r);
  f+=log_sigma+0.5*square(r/sigma);


SEPARABLE_FUNCTION void f2(int i, const prevariable& log_sigma, const prevariable& log_p1, const prevariable& log_p2, const prevariable& Ai, const prevariable& Bi, const prevariable& Ai1, const prevariable& Bi1, const prevariable& Phi1,const prevariable& Phi2)

  dvariable sigma=exp(log_sigma);
  dvariable p1=exp(log_p1);
  dvariable p2=exp(log_p2);

  dvariable r=Y(i)-Ai-Bi*Y(i-1);
  resids(i)=value(r);
  f+=log_sigma+0.5*square(r/sigma);
  dvariable nu_1i=(Ai-Phi1*Ai1);
  dvariable nu_2i=(Bi-Phi2*Bi1);
  f+=log_p1+0.5*square(nu_1i/p1);
  f+=log_p2+0.5*square(nu_2i/p2);

There is a "gotcha" with this approach.
Note that one can set A(i) such that
Y(i)-A(i)-B(i)*Y(i-1)=0

Then if you let sigma->0 the log-likelihood -> infinity.
So the "real" answer is a local maximum. You can stabilize the
estimation by putting a lower bound on sigma or keeping sigma fixed
at a reasonable value for a while. 
A: Would it possible to treat Y(t-1) as an exogenous variable, and estimate this state space model with kalman filter. Then the routine for estimating the coefficient is very standard.
