[pystatsmodels] Stochastic regression coeficients

[Alastair Heggie]

I would like to create a model in which the regression coefficients for
multiple exogenous covariates evolve with a seasonal component and a random
walk component.

For a single exogenous covariate I can create a local level model with a
seasonal component and modify the design matrix to use the exogenous
covariate like this:

dynamic_mod=sm.tsa.UnobservedComponents(endog,"local level",
freq_seasonal=
[{"period":24,"harmonics":nharmonics}])

exog_seasonal_array=np.vstack(tuple(np.array([exog,np.repeat(0,nobs)])
for i in range(0,nharmonics)))

dynamic_mod.ssm['design'] = np.asfortranarray(
np.array([np.vstack((exog,exog_seasonal_array))]))
dynamic_mod.ssm._time_invariant=False

I can use a similar approach to give the coefficients of multiple
covariates a seasonal component by just adding more components to
freq_seasonal. However, since the UnobservedComponents class only supports
a single local level component I can't use this method to have multiple
regression coefficients evolve according to a random walk.

In order to achieve what I want I think I will either need to create a new
class extending MLEModel or modify UnobservedComponents to allow stochastic
regression coefficients.

So far I have attempted the former: writing my own class that extends
MLEModel which I call SDR (Seasonal Dynamic Regerssion). However I am
having problems with convergence.

At the moment I set the starting parameters all as zeros, with the
exception of the variance of the irregular term which I set to the variance
of the data. However, I don't think the starting values are the root of the
problem. I can check if the starting values matter in the case with only 1
exogenous covariate because I can model this using the UnobservedComponents
class (as shown in the code block above). I use the optimised parameter
values for the starting parameters in my SDR class. When I do this the max
likelihood optimisation still fails to converge. However, the smoothed
state estimates from these starting parameters are the same as from the
UnobservedComponents model, which gives me more confidence that I have at
least defined the state space matrices correctly.

I wonder if part of the problem could be that I do not constrain the
parameters as in UnobservedComponents using a transform_params method.

I am hoping someone can help with these questions:

1) Is there is any fundamental problem with the kind of model I am trying
to form?
2) Is there is an easier way to implement this which I have not thought of
that does not involve writing a new class?
3) Are better starting parameters likely to help?
4) What might be responsible for the convergence problems?
5) Assuming I am able to implement this, would a pull request to add
stochastic regression components to the UnobservedComponents class be
welcomed?

I've attached my SDR class in case anyone is able to take a look.

Thanks!

[Chad Fulton]

Post by Alastair Heggie
I would like to create a model in which the regression coefficients for
multiple exogenous covariates evolve with a seasonal component and a random
walk component.
For a single exogenous covariate I can create a local level model with a
seasonal component and modify the design matrix to use the exogenous
dynamic_mod=sm.tsa.UnobservedComponents(endog,"local level",
freq_seasonal=
[{"period":24,"harmonics":nharmonics}])
exog_seasonal_array=np.vstack(tuple(np.array([exog,np.repeat(0,nobs)])
for i in range(0,nharmonics)))
dynamic_mod.ssm['design'] = np.asfortranarray(
np.array([np.vstack((exog,exog_seasonal_array))]))
dynamic_mod.ssm._time_invariant=False
I can use a similar approach to give the coefficients of multiple
covariates a seasonal component by just adding more components to
freq_seasonal. However, since the UnobservedComponents class only supports
a single local level component I can't use this method to have multiple
regression coefficients evolve according to a random walk.
In order to achieve what I want I think I will either need to create a new
class extending MLEModel or modify UnobservedComponents to allow stochastic
regression coefficients.
So far I have attempted the former: writing my own class that extends
MLEModel which I call SDR (Seasonal Dynamic Regerssion). However I am
having problems with convergence.
At the moment I set the starting parameters all as zeros, with the
exception of the variance of the irregular term which I set to the variance
of the data. However, I don't think the starting values are the root of the
problem. I can check if the starting values matter in the case with only 1
exogenous covariate because I can model this using the UnobservedComponents
class (as shown in the code block above). I use the optimised parameter
values for the starting parameters in my SDR class. When I do this the max
likelihood optimisation still fails to converge. However, the smoothed
state estimates from these starting parameters are the same as from the
UnobservedComponents model, which gives me more confidence that I have at
least defined the state space matrices correctly.
I wonder if part of the problem could be that I do not constrain the
parameters as in UnobservedComponents using a transform_params method.

Just to make sure I've understood the model:

y_t = beta_1,t x_1,t + beta_2,t x_2,t + eps_t

beta_1,t = mu_1,t + gamma_1,t
mu_1,t = mu_1,t-1 + eta_1,t

beta_2,t = mu_2,t + gamma_2,t + eta_2,t
mu_2,t = mu_2,t-1 + eta_1,t

and where gamma_1,t and gamma_2,t are separate, maybe stochastic, cycles.
Is that right?

*Best guess at the problem*

Regardless, I think the main problem is as you suspected, that you're not
using transform_params. If you add a `print(variances)` into your update
method and run `fit`, you can see that the optimizer quickly tries negative
values for the parameters, and the loglike function returns `nan`.

A quick fix is just to square the parameters in the update function (e.g.
add `params = params**2` right after the line `params = super(SDR,
self).update(params, **kwargs)`), but if you do this then you'll also need
to square the estimated parameters to get the estimates of the variance.
It's probably better in the long run to use transform_params.

*A secondary problem*

Something else I noticed is that the model you're setting up isn't like the
model I described above (although you can ignore this comment if that's not
the model you have in mind). For example, if I do:

mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4}])
mod['design', ..., 0]

then I get:

array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ]])

but for the model I wrote above, I think you would want (assuming you're
using the frequency seasonals):

array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ,
2710.349, 0. , 2710.349]])

*Your questions*

1) Is there is any fundamental problem with the kind of model I am trying

Post by Alastair Heggie
to form?

Nope, I think it should be fine.

Post by Alastair Heggie
2) Is there is an easier way to implement this which I have not thought of
that does not involve writing a new class?

I don't think so, we don't have anything for time-varying coefficients
other than simple random-walk coefficients in SARIMAX.

Post by Alastair Heggie
3) Are better starting parameters likely to help?

Always a good idea, but probably not the main problem here.

Post by Alastair Heggie
4) What might be responsible for the convergence problems?

See above.

Post by Alastair Heggie
5) Assuming I am able to implement this, would a pull request to add
stochastic regression components to the UnobservedComponents class be
welcomed?

Absolutely. How exactly we'd integrate it would depend on what you wanted
to actually put in the PR.

- If you want to stick within the UnobservedComponents class, then I think
we wouldn't want to go farther than the usual random walk coefficients. In
fact, we already allow the special case in which the random walk variance
is equal to zero (i.e. non-time-varying coefficients) with the argument
mle_regression=False, so extending that to allow time-varying coefficients
should be pretty straightforward (the main thing is adding the variance
parameters).

- If you are thinking about a PR with the full the case of regression with
time-varying coefficients that include both a random walk and seasonal
component, then I would think that might make more sense as a separate
class.

Best,
Chad

[Alastair Heggie]

Thanks for your comments on this Chad. Adding transform_params and
untransform_params methods to square the parameters works as desired and
I'm getting some nice results. I've attached the revised class for
interests sake.

Your understanding of the model is correct but I think my class sets up the
model correctly if you do:
mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4},{'period':4}])
(i.e. specifying the harmonics for each coefficient separately)

I'll have a think about pull requests to incorporate this into
statsmodels.I think it would take quite a bit of work to implement the full
the case of regression with time-varying coefficients but it would be an
interesting exercise for me. I suppose I would tackle this by trying to
give my revised SDR class a similar structure to Unobserved Components. Out
of interest, what do you think such as class should be called? Possibly
DynamicRegressionCoeficients? In any case, implementing random walk
coefficients in UnobservedComponents seems like a good place to start.

Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I would like to create a model in which the regression coefficients for
multiple exogenous covariates evolve with a seasonal component and a random
walk component.
For a single exogenous covariate I can create a local level model with a
seasonal component and modify the design matrix to use the exogenous
dynamic_mod=sm.tsa.UnobservedComponents(endog,"local level",
freq_seasonal=
[{"period":24,"harmonics":nharmonics}])
exog_seasonal_array=np.vstack(tuple(np.array([exog,np.repeat(0,nobs)])
for i in range(0,nharmonics)))
dynamic_mod.ssm['design'] = np.asfortranarray(
np.array([np.vstack((exog,exog_seasonal_array))]))
dynamic_mod.ssm._time_invariant=False
I can use a similar approach to give the coefficients of multiple
covariates a seasonal component by just adding more components to
freq_seasonal. However, since the UnobservedComponents class only supports
a single local level component I can't use this method to have multiple
regression coefficients evolve according to a random walk.
In order to achieve what I want I think I will either need to create a
new class extending MLEModel or modify UnobservedComponents to allow
stochastic regression coefficients.
So far I have attempted the former: writing my own class that extends
MLEModel which I call SDR (Seasonal Dynamic Regerssion). However I am
having problems with convergence.
At the moment I set the starting parameters all as zeros, with the
exception of the variance of the irregular term which I set to the variance
of the data. However, I don't think the starting values are the root of the
problem. I can check if the starting values matter in the case with only 1
exogenous covariate because I can model this using the UnobservedComponents
class (as shown in the code block above). I use the optimised parameter
values for the starting parameters in my SDR class. When I do this the max
likelihood optimisation still fails to converge. However, the smoothed
state estimates from these starting parameters are the same as from the
UnobservedComponents model, which gives me more confidence that I have at
least defined the state space matrices correctly.
I wonder if part of the problem could be that I do not constrain the
parameters as in UnobservedComponents using a transform_params method.

y_t = beta_1,t x_1,t + beta_2,t x_2,t + eps_t
beta_1,t = mu_1,t + gamma_1,t
mu_1,t = mu_1,t-1 + eta_1,t
beta_2,t = mu_2,t + gamma_2,t + eta_2,t
mu_2,t = mu_2,t-1 + eta_1,t
and where gamma_1,t and gamma_2,t are separate, maybe stochastic, cycles.
Is that right?
*Best guess at the problem*
Regardless, I think the main problem is as you suspected, that you're not
using transform_params. If you add a `print(variances)` into your update
method and run `fit`, you can see that the optimizer quickly tries negative
values for the parameters, and the loglike function returns `nan`.
A quick fix is just to square the parameters in the update function (e.g.
add `params = params**2` right after the line `params = super(SDR,
self).update(params, **kwargs)`), but if you do this then you'll also need
to square the estimated parameters to get the estimates of the variance.
It's probably better in the long run to use transform_params.
*A secondary problem*
Something else I noticed is that the model you're setting up isn't like
the model I described above (although you can ignore this comment if that's
mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4}])
mod['design', ..., 0]
array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ]])
but for the model I wrote above, I think you would want (assuming you're
array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ,
2710.349, 0. , 2710.349]])
*Your questions*
1) Is there is any fundamental problem with the kind of model I am trying

Post by Alastair Heggie
to form?

Nope, I think it should be fine.

Post by Alastair Heggie
2) Is there is an easier way to implement this which I have not thought
of that does not involve writing a new class?

I don't think so, we don't have anything for time-varying coefficients
other than simple random-walk coefficients in SARIMAX.

Post by Alastair Heggie
3) Are better starting parameters likely to help?

Always a good idea, but probably not the main problem here.

Post by Alastair Heggie
4) What might be responsible for the convergence problems?

See above.

Post by Alastair Heggie
5) Assuming I am able to implement this, would a pull request to add
stochastic regression components to the UnobservedComponents class be
welcomed?

Absolutely. How exactly we'd integrate it would depend on what you wanted
to actually put in the PR.
- If you want to stick within the UnobservedComponents class, then I think
we wouldn't want to go farther than the usual random walk coefficients. In
fact, we already allow the special case in which the random walk variance
is equal to zero (i.e. non-time-varying coefficients) with the argument
mle_regression=False, so extending that to allow time-varying coefficients
should be pretty straightforward (the main thing is adding the variance
parameters).
- If you are thinking about a PR with the full the case of regression with
time-varying coefficients that include both a random walk and seasonal
component, then I would think that might make more sense as a separate
class.
Best,
Chad

[Alastair Heggie]

I'd like to return to this topic and explore what would be necessary to
integrate my work on stochastic regression coefficients into Statsmodels.

The following gist contains a new class, SDR.py and an ipynb demonstrating
my work on some simulated data:
https://gist.github.com/AliHeggie/c6ef14e5f03f57218484888b3a2ceee8

The class SDR.py implements the model discussed above where the regression
coefficients evolve as a random walk plus stochastic trigonometric cycles.

I was hoping someone could help me assemble a To Do list of style changes
and feature additions that be necessary to get this into Statsmodels.

Some things that spring to mind:
- Tests
- A new results class with plotting methods (At the moment I have a plot
method in SDR.py but it is not associated with any class, it just takes a
MLEResultsWrapper as an argument)
- I would like to add the option to specify and ARMA model for the
observation errors.
- More possibilities for the model of the regression coefficients, dummy
seasonal variables instead of trigonometric for example.

Best wishes,
Alastair

Post by Alastair Heggie
Thanks for your comments on this Chad. Adding transform_params and
untransform_params methods to square the parameters works as desired and
I'm getting some nice results. I've attached the revised class for
interests sake.
Your understanding of the model is correct but I think my class sets up
mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4},{'period':4
}])
(i.e. specifying the harmonics for each coefficient separately)
I'll have a think about pull requests to incorporate this into
statsmodels.I think it would take quite a bit of work to implement the full
the case of regression with time-varying coefficients but it would be an
interesting exercise for me. I suppose I would tackle this by trying to
give my revised SDR class a similar structure to Unobserved Components. Out
of interest, what do you think such as class should be called? Possibly
DynamicRegressionCoeficients? In any case, implementing random walk
coefficients in UnobservedComponents seems like a good place to start.
Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I would like to create a model in which the regression coefficients for
multiple exogenous covariates evolve with a seasonal component and a random
walk component.
For a single exogenous covariate I can create a local level model with a
seasonal component and modify the design matrix to use the exogenous
dynamic_mod=sm.tsa.UnobservedComponents(endog,"local level",
freq_seasonal=
[{"period":24,"harmonics":nharmonics}])
exog_seasonal_array=np.vstack(tuple(np.array([exog,np.repeat(0,nobs)])
for i in range(0,nharmonics)))
dynamic_mod.ssm['design'] = np.asfortranarray(
np.array([np.vstack((exog,exog_seasonal_array))]))
dynamic_mod.ssm._time_invariant=False
I can use a similar approach to give the coefficients of multiple
covariates a seasonal component by just adding more components to
freq_seasonal. However, since the UnobservedComponents class only supports
a single local level component I can't use this method to have multiple
regression coefficients evolve according to a random walk.
In order to achieve what I want I think I will either need to create a
new class extending MLEModel or modify UnobservedComponents to allow
stochastic regression coefficients.
So far I have attempted the former: writing my own class that extends
MLEModel which I call SDR (Seasonal Dynamic Regerssion). However I am
having problems with convergence.
At the moment I set the starting parameters all as zeros, with the
exception of the variance of the irregular term which I set to the variance
of the data. However, I don't think the starting values are the root of the
problem. I can check if the starting values matter in the case with only 1
exogenous covariate because I can model this using the UnobservedComponents
class (as shown in the code block above). I use the optimised parameter
values for the starting parameters in my SDR class. When I do this the max
likelihood optimisation still fails to converge. However, the smoothed
state estimates from these starting parameters are the same as from the
UnobservedComponents model, which gives me more confidence that I have at
least defined the state space matrices correctly.
I wonder if part of the problem could be that I do not constrain the
parameters as in UnobservedComponents using a transform_params method.

y_t = beta_1,t x_1,t + beta_2,t x_2,t + eps_t
beta_1,t = mu_1,t + gamma_1,t
mu_1,t = mu_1,t-1 + eta_1,t
beta_2,t = mu_2,t + gamma_2,t + eta_2,t
mu_2,t = mu_2,t-1 + eta_1,t
and where gamma_1,t and gamma_2,t are separate, maybe stochastic, cycles.
Is that right?
*Best guess at the problem*
Regardless, I think the main problem is as you suspected, that you're not
using transform_params. If you add a `print(variances)` into your update
method and run `fit`, you can see that the optimizer quickly tries negative
values for the parameters, and the loglike function returns `nan`.
A quick fix is just to square the parameters in the update function (e.g.
add `params = params**2` right after the line `params = super(SDR,
self).update(params, **kwargs)`), but if you do this then you'll also need
to square the estimated parameters to get the estimates of the variance.
It's probably better in the long run to use transform_params.
*A secondary problem*
Something else I noticed is that the model you're setting up isn't like
the model I described above (although you can ignore this comment if that's
mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4}])
mod['design', ..., 0]
array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ]])
but for the model I wrote above, I think you would want (assuming you're
array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ,
2710.349, 0. , 2710.349]])
*Your questions*
1) Is there is any fundamental problem with the kind of model I am trying

Post by Alastair Heggie
to form?

Nope, I think it should be fine.

Post by Alastair Heggie
2) Is there is an easier way to implement this which I have not thought
of that does not involve writing a new class?

I don't think so, we don't have anything for time-varying coefficients
other than simple random-walk coefficients in SARIMAX.

Post by Alastair Heggie
3) Are better starting parameters likely to help?

Always a good idea, but probably not the main problem here.

Post by Alastair Heggie
4) What might be responsible for the convergence problems?

See above.

Post by Alastair Heggie
5) Assuming I am able to implement this, would a pull request to add
stochastic regression components to the UnobservedComponents class be
welcomed?

Absolutely. How exactly we'd integrate it would depend on what you wanted
to actually put in the PR.
- If you want to stick within the UnobservedComponents class, then I
think we wouldn't want to go farther than the usual random walk
coefficients. In fact, we already allow the special case in which the
random walk variance is equal to zero (i.e. non-time-varying coefficients)
with the argument mle_regression=False, so extending that to allow
time-varying coefficients should be pretty straightforward (the main thing
is adding the variance parameters).
- If you are thinking about a PR with the full the case of regression
with time-varying coefficients that include both a random walk and seasonal
component, then I would think that might make more sense as a separate
class.
Best,
Chad

[Alastair Heggie]

As a follow up to this I have an alternative proposition for which I'd be
interested to get feedback. I think this would be a much more general
approach that I think should be relatively easy to implement using existing
classes.

For any state space model that allows exogenous covariates, we could allow
as an additional kwarg a list of state space models which define how the
coefficients of each of the exogenous covariates evolves. We can then
append the system matrices of the model for the exogenous coefficient to
the system matrices for the top level model to create a new combined model.

For example, suppose we want to define a model with ARIMA errors and a
single exogenous regressor who's coefficient evolves according to a
structural model. We can use the SARIMAX model as our base class and supply
the structural model as an argument.

- We would create the design matrix of the combined model by appending to
the bottom of the SARIMAX design matrix the design matrix of the structural
model multiplied by the exogenous regressor at every time step.
- We would create the transition matrix of the combined model by combining
the transition matrices of the two models as a block diagonal matrix.
- Similarly the state covariance matrix is a block diagonal matrix created
from the SARIMAX and structural state covariance matrix. This means that to
update the model with new parameters we can call the structural model's
update method and then use the state covariance matrix of that updated
model to rebuild the state_covariance model of the new model.

I've had a crack at modifying the sarimax class to implementing this here:
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegression.py.
At the moment something goes wrong when we try to fit the resulting model
as you can see in the last cell of this notebook:
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegressionTests.ipynb

I'd be interested to hear if anyone thinks this approach is worth pursuing
further to get a fairly general class of dynamic regression model.

Best wishes,
Alastair

Post by Alastair Heggie
I'd like to return to this topic and explore what would be necessary to
integrate my work on stochastic regression coefficients into Statsmodels.
The following gist contains a new class, SDR.py and an ipynb demonstrating
https://gist.github.com/AliHeggie/c6ef14e5f03f57218484888b3a2ceee8
The class SDR.py implements the model discussed above where the regression
coefficients evolve as a random walk plus stochastic trigonometric cycles.
I was hoping someone could help me assemble a To Do list of style changes
and feature additions that be necessary to get this into Statsmodels.
- Tests
- A new results class with plotting methods (At the moment I have a plot
method in SDR.py but it is not associated with any class, it just takes a
MLEResultsWrapper as an argument)
- I would like to add the option to specify and ARMA model for the
observation errors.
- More possibilities for the model of the regression coefficients, dummy
seasonal variables instead of trigonometric for example.
Best wishes,
Alastair

Post by Alastair Heggie
Thanks for your comments on this Chad. Adding transform_params and
untransform_params methods to square the parameters works as desired and
I'm getting some nice results. I've attached the revised class for
interests sake.
Your understanding of the model is correct but I think my class sets up
mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4},{'period':4
}])
(i.e. specifying the harmonics for each coefficient separately)
I'll have a think about pull requests to incorporate this into
statsmodels.I think it would take quite a bit of work to implement the full
the case of regression with time-varying coefficients but it would be an
interesting exercise for me. I suppose I would tackle this by trying to
give my revised SDR class a similar structure to Unobserved Components. Out
of interest, what do you think such as class should be called? Possibly
DynamicRegressionCoeficients? In any case, implementing random walk
coefficients in UnobservedComponents seems like a good place to start.
Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
I would like to create a model in which the regression coefficients for
multiple exogenous covariates evolve with a seasonal component and a random
walk component.
For a single exogenous covariate I can create a local level model with
a seasonal component and modify the design matrix to use the exogenous
dynamic_mod=sm.tsa.UnobservedComponents(endog,"local level",
freq_seasonal=
[{"period":24,"harmonics":nharmonics}])
exog_seasonal_array=np.vstack(tuple(np.array([exog,np.repeat(0,nobs)])
for i in range(0,nharmonics)))
dynamic_mod.ssm['design'] = np.asfortranarray(
np.array([np.vstack((exog,exog_seasonal_array))]))
dynamic_mod.ssm._time_invariant=False
I can use a similar approach to give the coefficients of multiple
covariates a seasonal component by just adding more components to
freq_seasonal. However, since the UnobservedComponents class only supports
a single local level component I can't use this method to have multiple
regression coefficients evolve according to a random walk.
In order to achieve what I want I think I will either need to create a
new class extending MLEModel or modify UnobservedComponents to allow
stochastic regression coefficients.
So far I have attempted the former: writing my own class that extends
MLEModel which I call SDR (Seasonal Dynamic Regerssion). However I am
having problems with convergence.
At the moment I set the starting parameters all as zeros, with the
exception of the variance of the irregular term which I set to the variance
of the data. However, I don't think the starting values are the root of the
problem. I can check if the starting values matter in the case with only 1
exogenous covariate because I can model this using the UnobservedComponents
class (as shown in the code block above). I use the optimised parameter
values for the starting parameters in my SDR class. When I do this the max
likelihood optimisation still fails to converge. However, the smoothed
state estimates from these starting parameters are the same as from the
UnobservedComponents model, which gives me more confidence that I have at
least defined the state space matrices correctly.
I wonder if part of the problem could be that I do not constrain the
parameters as in UnobservedComponents using a transform_params method.

y_t = beta_1,t x_1,t + beta_2,t x_2,t + eps_t
beta_1,t = mu_1,t + gamma_1,t
mu_1,t = mu_1,t-1 + eta_1,t
beta_2,t = mu_2,t + gamma_2,t + eta_2,t
mu_2,t = mu_2,t-1 + eta_1,t
and where gamma_1,t and gamma_2,t are separate, maybe stochastic,
cycles. Is that right?
*Best guess at the problem*
Regardless, I think the main problem is as you suspected, that you're
not using transform_params. If you add a `print(variances)` into your
update method and run `fit`, you can see that the optimizer quickly tries
negative values for the parameters, and the loglike function returns `nan`.
A quick fix is just to square the parameters in the update function
(e.g. add `params = params**2` right after the line `params = super(SDR,
self).update(params, **kwargs)`), but if you do this then you'll also need
to square the estimated parameters to get the estimates of the variance.
It's probably better in the long run to use transform_params.
*A secondary problem*
Something else I noticed is that the model you're setting up isn't like
the model I described above (although you can ignore this comment if that's
mod = SDR(dta['infl'], dta[['cpi', 'realgdp']], [{'period':4}])
mod['design', ..., 0]
array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ]])
but for the model I wrote above, I think you would want (assuming you're
array([[ 28.98 , 2710.349, 28.98 , 0. , 28.98 , 0. ,
2710.349, 0. , 2710.349]])
*Your questions*
1) Is there is any fundamental problem with the kind of model I am

Post by Alastair Heggie
trying to form?

Nope, I think it should be fine.

Post by Alastair Heggie
2) Is there is an easier way to implement this which I have not thought
of that does not involve writing a new class?

I don't think so, we don't have anything for time-varying coefficients
other than simple random-walk coefficients in SARIMAX.

Post by Alastair Heggie
3) Are better starting parameters likely to help?

Always a good idea, but probably not the main problem here.

Post by Alastair Heggie
4) What might be responsible for the convergence problems?

See above.

Post by Alastair Heggie
5) Assuming I am able to implement this, would a pull request to add
stochastic regression components to the UnobservedComponents class be
welcomed?

Absolutely. How exactly we'd integrate it would depend on what you
wanted to actually put in the PR.
- If you want to stick within the UnobservedComponents class, then I
think we wouldn't want to go farther than the usual random walk
coefficients. In fact, we already allow the special case in which the
random walk variance is equal to zero (i.e. non-time-varying coefficients)
with the argument mle_regression=False, so extending that to allow
time-varying coefficients should be pretty straightforward (the main thing
is adding the variance parameters).
- If you are thinking about a PR with the full the case of regression
with time-varying coefficients that include both a random walk and seasonal
component, then I would think that might make more sense as a separate
class.
Best,
Chad

[Chad Fulton]

Post by Alastair Heggie
As a follow up to this I have an alternative proposition for which I'd be
interested to get feedback. I think this would be a much more general
approach that I think should be relatively easy to implement using existing
classes.
For any state space model that allows exogenous covariates, we could allow
as an additional kwarg a list of state space models which define how the
coefficients of each of the exogenous covariates evolves. We can then
append the system matrices of the model for the exogenous coefficient to
the system matrices for the top level model to create a new combined model.
For example, suppose we want to define a model with ARIMA errors and a
single exogenous regressor who's coefficient evolves according to a
structural model. We can use the SARIMAX model as our base class and supply
the structural model as an argument.
- We would create the design matrix of the combined model by appending to
the bottom of the SARIMAX design matrix the design matrix of the structural
model multiplied by the exogenous regressor at every time step.
- We would create the transition matrix of the combined model by combining
the transition matrices of the two models as a block diagonal matrix.
- Similarly the state covariance matrix is a block diagonal matrix created
from the SARIMAX and structural state covariance matrix. This means that to
update the model with new parameters we can call the structural model's
update method and then use the state covariance matrix of that updated
model to rebuild the state_covariance model of the new model.
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegression.py.
At the moment something goes wrong when we try to fit the resulting model
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegressionTests.ipynb
I'd be interested to hear if anyone thinks this approach is worth pursuing
further to get a fairly general class of dynamic regression model.
Best wishes,
Alastair

Post by Alastair Heggie
I'd like to return to this topic and explore what would be necessary to
integrate my work on stochastic regression coefficients into Statsmodels.
The following gist contains a new class, SDR.py and an ipynb
https://gist.github.com/AliHeggie/c6ef14e5f03f57218484888b3a2ceee8
The class SDR.py implements the model discussed above where the
regression coefficients evolve as a random walk plus stochastic
trigonometric cycles.
I was hoping someone could help me assemble a To Do list of style changes
and feature additions that be necessary to get this into Statsmodels.
- Tests
- A new results class with plotting methods (At the moment I have a plot
method in SDR.py but it is not associated with any class, it just takes a
MLEResultsWrapper as an argument)
- I would like to add the option to specify and ARMA model for the
observation errors.
- More possibilities for the model of the regression coefficients, dummy
seasonal variables instead of trigonometric for example.
Best wishes,
Alastair

Hi Alastair,

Thank you for sharing all of your work on this! It looks very promising,
and I'll look it over your SDR class in detail as soon as I can. I'd be
happy to provide a TODO list with some details on how we could get it in
shape for being merged.

One thing that would be really interesting to me (will not be necessary for
merge, but would be nice) - is there a public dataset which can demonstrate
the use of this model? If so, it would be really neat to see it in action,
if you could post a short tutorial notebook on actual data.

-----------

For your more recent e-mail, that is a very intriguing idea for building
more complex state space models. I'll take a closer look at what you've
done soon, and see if I can see what is left to do to make it fully
workable.

Thanks again for posting all of this!
Chad

[Alastair Heggie]

Hi Chad,

I've just pushed some changes to the SDR class in my DynamicRegression repo
to implement ARMA errors and non dynamic regression coefficients (i.e. 0
covariance, but still potentially seasonal).

I would like to make public the data from the project that motivated this
work. The data is from industry so I need to OK everything with my partners
first however.

Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
As a follow up to this I have an alternative proposition for which I'd be
interested to get feedback. I think this would be a much more general
approach that I think should be relatively easy to implement using existing
classes.
For any state space model that allows exogenous covariates, we could
allow as an additional kwarg a list of state space models which define how
the coefficients of each of the exogenous covariates evolves. We can then
append the system matrices of the model for the exogenous coefficient to
the system matrices for the top level model to create a new combined model.
For example, suppose we want to define a model with ARIMA errors and a
single exogenous regressor who's coefficient evolves according to a
structural model. We can use the SARIMAX model as our base class and supply
the structural model as an argument.
- We would create the design matrix of the combined model by appending to
the bottom of the SARIMAX design matrix the design matrix of the structural
model multiplied by the exogenous regressor at every time step.
- We would create the transition matrix of the combined model by
combining the transition matrices of the two models as a block diagonal
matrix.
- Similarly the state covariance matrix is a block diagonal matrix
created from the SARIMAX and structural state covariance matrix. This means
that to update the model with new parameters we can call the structural
model's update method and then use the state covariance matrix of that
updated model to rebuild the state_covariance model of the new model.
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegression.py.
At the moment something goes wrong when we try to fit the resulting model
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegressionTests.ipynb
I'd be interested to hear if anyone thinks this approach is worth
pursuing further to get a fairly general class of dynamic regression model.
Best wishes,
Alastair

Post by Alastair Heggie
I'd like to return to this topic and explore what would be necessary to
integrate my work on stochastic regression coefficients into Statsmodels.
The following gist contains a new class, SDR.py and an ipynb
https://gist.github.com/AliHeggie/c6ef14e5f03f57218484888b3a2ceee8
The class SDR.py implements the model discussed above where the
regression coefficients evolve as a random walk plus stochastic
trigonometric cycles.
I was hoping someone could help me assemble a To Do list of style
changes and feature additions that be necessary to get this into
Statsmodels.
- Tests
- A new results class with plotting methods (At the moment I have a plot
method in SDR.py but it is not associated with any class, it just takes a
MLEResultsWrapper as an argument)
- I would like to add the option to specify and ARMA model for the
observation errors.
- More possibilities for the model of the regression coefficients, dummy
seasonal variables instead of trigonometric for example.
Best wishes,
Alastair

Hi Alastair,
Thank you for sharing all of your work on this! It looks very promising,
and I'll look it over your SDR class in detail as soon as I can. I'd be
happy to provide a TODO list with some details on how we could get it in
shape for being merged.
One thing that would be really interesting to me (will not be necessary
for merge, but would be nice) - is there a public dataset which can
demonstrate the use of this model? If so, it would be really neat to see it
in action, if you could post a short tutorial notebook on actual data.
-----------
For your more recent e-mail, that is a very intriguing idea for building
more complex state space models. I'll take a closer look at what you've
done soon, and see if I can see what is left to do to make it fully
workable.
Thanks again for posting all of this!
Chad

[Alastair Heggie]

I've forked the statsmodels repo and moved the development of this feature
to a new branch called dynamic-regression:
https://github.com/AliHeggie/statsmodels/tree/dynamic-regression.

Using http://www.statsmodels.org/stable/dev/ as a guide to what do do to
get the code in a fit state for merge.

One thing I'm not quite clear on is what is the appropriate stage in
development to submit a pull request. Should I wait till I think the code
is totally up to standard or just use a PR as a way to invite feedback on
code that is probably not yet ready for merging?

Post by Alastair Heggie
Hi Chad,
I've just pushed some changes to the SDR class in my DynamicRegression
repo to implement ARMA errors and non dynamic regression coefficients (i.e.
0 covariance, but still potentially seasonal).
I would like to make public the data from the project that motivated this
work. The data is from industry so I need to OK everything with my partners
first however.
Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
As a follow up to this I have an alternative proposition for which I'd
be interested to get feedback. I think this would be a much more general
approach that I think should be relatively easy to implement using existing
classes.
For any state space model that allows exogenous covariates, we could
allow as an additional kwarg a list of state space models which define how
the coefficients of each of the exogenous covariates evolves. We can then
append the system matrices of the model for the exogenous coefficient to
the system matrices for the top level model to create a new combined model.
For example, suppose we want to define a model with ARIMA errors and a
single exogenous regressor who's coefficient evolves according to a
structural model. We can use the SARIMAX model as our base class and supply
the structural model as an argument.
- We would create the design matrix of the combined model by appending
to the bottom of the SARIMAX design matrix the design matrix of the
structural model multiplied by the exogenous regressor at every time step.
- We would create the transition matrix of the combined model by
combining the transition matrices of the two models as a block diagonal
matrix.
- Similarly the state covariance matrix is a block diagonal matrix
created from the SARIMAX and structural state covariance matrix. This means
that to update the model with new parameters we can call the structural
model's update method and then use the state covariance matrix of that
updated model to rebuild the state_covariance model of the new model.
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegression.py.
At the moment something goes wrong when we try to fit the resulting model
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegressionTests.ipynb
I'd be interested to hear if anyone thinks this approach is worth
pursuing further to get a fairly general class of dynamic regression model.
Best wishes,
Alastair

Post by Alastair Heggie
I'd like to return to this topic and explore what would be necessary to
integrate my work on stochastic regression coefficients into Statsmodels.
The following gist contains a new class, SDR.py and an ipynb
https://gist.github.com/AliHeggie/c6ef14e5f03f57218484888b3a2ceee8
The class SDR.py implements the model discussed above where the
regression coefficients evolve as a random walk plus stochastic
trigonometric cycles.
I was hoping someone could help me assemble a To Do list of style
changes and feature additions that be necessary to get this into
Statsmodels.
- Tests
- A new results class with plotting methods (At the moment I have a
plot method in SDR.py but it is not associated with any class, it just
takes a MLEResultsWrapper as an argument)
- I would like to add the option to specify and ARMA model for the
observation errors.
- More possibilities for the model of the regression coefficients,
dummy seasonal variables instead of trigonometric for example.
Best wishes,
Alastair

Hi Alastair,
Thank you for sharing all of your work on this! It looks very promising,
and I'll look it over your SDR class in detail as soon as I can. I'd be
happy to provide a TODO list with some details on how we could get it in
shape for being merged.
One thing that would be really interesting to me (will not be necessary
for merge, but would be nice) - is there a public dataset which can
demonstrate the use of this model? If so, it would be really neat to see it
in action, if you could post a short tutorial notebook on actual data.
-----------
For your more recent e-mail, that is a very intriguing idea for building
more complex state space models. I'll take a closer look at what you've
done soon, and see if I can see what is left to do to make it fully
workable.
Thanks again for posting all of this!
Chad

[Chad Fulton]

Post by Alastair Heggie
I've forked the statsmodels repo and moved the development of this feature
https://github.com/AliHeggie/statsmodels/tree/dynamic-regression.
Using http://www.statsmodels.org/stable/dev/ as a guide to what do do to
get the code in a fit state for merge.
One thing I'm not quite clear on is what is the appropriate stage in
development to submit a pull request. Should I wait till I think the code
is totally up to standard or just use a PR as a way to invite feedback on
code that is probably not yet ready for merging?

Hi Alastair,

Apologies for the delay. I've now had a chance to take a first pass at your
SDR class and I think it looks quite nice. Some notes / questions:

- Please do submit a PR with the code whenever convenient, that is the best
way for us to e.g. add comments.

- It would be helpful if you could document the model that you have in mind
in the class docstring. If you have a reference for this model (e.g.
showing that this model is important / useful in the literature), that
would be the best case. (Of course the basic random walk time-varying
parameters regression model is well used, but I wonder if there is a
particular reference for the random walk + seasonal harmonics + ARMA form
of the parameter transition).

- It would be helpful if you document the state space form of the model
(for example, I see that you're using what we've called "Hamilton's" form
for the ARMA component)

There are two things that come to mind that we should decide:

1. We need to decide "where to put" / "how to describe" this class. A basic
random walk time-varying parameters regression would definitely fit in as a
primary `tsa.statespace` class, but we haven't yet thought about what to do
with more ad-hoc model definitions. One option might be to split this into
two classes - one with the basic TVP model and one extending it that adds
support for the additional types of coefficient dynamics.

2. This model is connected with our `RecursiveLS` class, since the basic
TVP model is essentially `RecursiveLS` with non-zero state error variances.
Should this class (or the base TVP class, if we end up splitting it) extend
the `RecursiveLS` class? (It is entirely possible that it should not extend
it, since the usual uses of the two models seems very different, and I'm
not sure how many of the `RecursiveLS` results, like the recursive
residuals, transfer over usefully).

Best,
Chad

[Alastair Heggie]

Thanks Chad,

I've submitted a pull request for a simplified version of the model that
just implements time varying coefficients that evolve according to a
local_level + trend + error UC model.

I'd like to just get this version up to standard to be included.

The question of where to put this is an interesting one. The current pull
request is for a basic random walk time-varying parameters regression so
would probably fit as a separate class. The extension I propose for this is
just to add more of the standard features of an unobserved components
models as optional features of coefficient transition dynamics. This seems
like a fairly natural extension to me and so makes sense as part of a basic
dynamic regression class, what do you think?

I'm not aware that the local level + seasonal harmonics + ARMA is seen in
the literature but I've found this version helpful in my research. When
it's ready I could provide a link to a working paper.

Post by Chad Fulton

Post by Alastair Heggie
I've forked the statsmodels repo and moved the development of this
https://github.com/AliHeggie/statsmodels/tree/dynamic-regression.
Using http://www.statsmodels.org/stable/dev/ as a guide to what do do to
get the code in a fit state for merge.
One thing I'm not quite clear on is what is the appropriate stage in
development to submit a pull request. Should I wait till I think the code
is totally up to standard or just use a PR as a way to invite feedback on
code that is probably not yet ready for merging?

Hi Alastair,
Apologies for the delay. I've now had a chance to take a first pass at
- Please do submit a PR with the code whenever convenient, that is the
best way for us to e.g. add comments.
- It would be helpful if you could document the model that you have in
mind in the class docstring. If you have a reference for this model (e.g.
showing that this model is important / useful in the literature), that
would be the best case. (Of course the basic random walk time-varying
parameters regression model is well used, but I wonder if there is a
particular reference for the random walk + seasonal harmonics + ARMA form
of the parameter transition).
- It would be helpful if you document the state space form of the model
(for example, I see that you're using what we've called "Hamilton's" form
for the ARMA component)
1. We need to decide "where to put" / "how to describe" this class. A
basic random walk time-varying parameters regression would definitely fit
in as a primary `tsa.statespace` class, but we haven't yet thought about
what to do with more ad-hoc model definitions. One option might be to split
this into two classes - one with the basic TVP model and one extending it
that adds support for the additional types of coefficient dynamics.
2. This model is connected with our `RecursiveLS` class, since the basic
TVP model is essentially `RecursiveLS` with non-zero state error variances.
Should this class (or the base TVP class, if we end up splitting it) extend
the `RecursiveLS` class? (It is entirely possible that it should not extend
it, since the usual uses of the two models seems very different, and I'm
not sure how many of the `RecursiveLS` results, like the recursive
residuals, transfer over usefully).
Best,
Chad

[Chad Fulton]

Post by Alastair Heggie
As a follow up to this I have an alternative proposition for which I'd be
interested to get feedback. I think this would be a much more general
approach that I think should be relatively easy to implement using existing
classes.
For any state space model that allows exogenous covariates, we could allow
as an additional kwarg a list of state space models which define how the
coefficients of each of the exogenous covariates evolves. We can then
append the system matrices of the model for the exogenous coefficient to
the system matrices for the top level model to create a new combined model.
For example, suppose we want to define a model with ARIMA errors and a
single exogenous regressor who's coefficient evolves according to a
structural model. We can use the SARIMAX model as our base class and supply
the structural model as an argument.
- We would create the design matrix of the combined model by appending to
the bottom of the SARIMAX design matrix the design matrix of the structural
model multiplied by the exogenous regressor at every time step.
- We would create the transition matrix of the combined model by combining
the transition matrices of the two models as a block diagonal matrix.
- Similarly the state covariance matrix is a block diagonal matrix created
from the SARIMAX and structural state covariance matrix. This means that to
update the model with new parameters we can call the structural model's
update method and then use the state covariance matrix of that updated
model to rebuild the state_covariance model of the new model.
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegression.py.
At the moment something goes wrong when we try to fit the resulting model
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegressionTests.ipynb
I'd be interested to hear if anyone thinks this approach is worth pursuing
further to get a fairly general class of dynamic regression model.
Best wishes,
Alastair

(I split this off from the other threat to keep things clear)

In general, I'm in favor of making our state space models more flexible,
but I have to say I'm not sure about this particular proposal, mainly for
two reasons:

1. It makes the already-complex SARIMAX class even more complex. In
addition to potentially adding bugs, there's a concern that it adds
difficult-to-debug corner cases, and it certainly adds to the already
existing problem that we don't test all possible combinations of each model.

2. Since state space models are so general, there are no end to the
possible additions we could make to each model, and this makes me think
that we should err towards the side of restricting our core state space
models to those that are "well-established" (i.e. models that some
reasonable record in the literature).

So two things come to mind:

a. We may want to establish a "contributed' section of state space models
that include models that are relatively less well-known.

b. I think a better long-term approach is to make state space models
compossible, something like the following:

mod_a = sm.tsa.SARIMAX(endog, order=(1, 0, 1))
mod_b = sm.tsa.UnobservedComponents(endog, 'llevel')

mod = mod_a + mod_b
res = mod.fit()

I've started laying some basic groundwork for this approach when I reworked
the state space initialization framework (
https://github.com/statsmodels/statsmodels/pull/4418). I didn't go so far
as to override the addition operator as in my pseudo-code here, but it does
allow for arbitrary nesting of initializations.

I'd be interested in any thoughts on this topic.

Best,
Chad

[Alastair Heggie]

I think I agree with all of the above. For the time-being I will focus my
efforts on a pull request for my dynamic regression class.

In the long term making state space models composable seems like a
desirable feature.

Best,
Alastair

Post by Chad Fulton

Post by Alastair Heggie
As a follow up to this I have an alternative proposition for which I'd be
interested to get feedback. I think this would be a much more general
approach that I think should be relatively easy to implement using existing
classes.
For any state space model that allows exogenous covariates, we could
allow as an additional kwarg a list of state space models which define how
the coefficients of each of the exogenous covariates evolves. We can then
append the system matrices of the model for the exogenous coefficient to
the system matrices for the top level model to create a new combined model.
For example, suppose we want to define a model with ARIMA errors and a
single exogenous regressor who's coefficient evolves according to a
structural model. We can use the SARIMAX model as our base class and supply
the structural model as an argument.
- We would create the design matrix of the combined model by appending to
the bottom of the SARIMAX design matrix the design matrix of the structural
model multiplied by the exogenous regressor at every time step.
- We would create the transition matrix of the combined model by
combining the transition matrices of the two models as a block diagonal
matrix.
- Similarly the state covariance matrix is a block diagonal matrix
created from the SARIMAX and structural state covariance matrix. This means
that to update the model with new parameters we can call the structural
model's update method and then use the state covariance matrix of that
updated model to rebuild the state_covariance model of the new model.
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegression.py.
At the moment something goes wrong when we try to fit the resulting model
https://github.com/AliHeggie/DynamicRegression/blob/master/DynamicRegressionTests.ipynb
I'd be interested to hear if anyone thinks this approach is worth
pursuing further to get a fairly general class of dynamic regression model.
Best wishes,
Alastair

(I split this off from the other threat to keep things clear)
In general, I'm in favor of making our state space models more flexible,
but I have to say I'm not sure about this particular proposal, mainly for
1. It makes the already-complex SARIMAX class even more complex. In
addition to potentially adding bugs, there's a concern that it adds
difficult-to-debug corner cases, and it certainly adds to the already
existing problem that we don't test all possible combinations of each model.
2. Since state space models are so general, there are no end to the
possible additions we could make to each model, and this makes me think
that we should err towards the side of restricting our core state space
models to those that are "well-established" (i.e. models that some
reasonable record in the literature).
a. We may want to establish a "contributed' section of state space models
that include models that are relatively less well-known.
b. I think a better long-term approach is to make state space models
mod_a = sm.tsa.SARIMAX(endog, order=(1, 0, 1))
mod_b = sm.tsa.UnobservedComponents(endog, 'llevel')
mod = mod_a + mod_b
res = mod.fit()
I've started laying some basic groundwork for this approach when I
reworked the state space initialization framework (
https://github.com/statsmodels/statsmodels/pull/4418). I didn't go so far
as to override the addition operator as in my pseudo-code here, but it does
allow for arbitrary nesting of initializations.
I'd be interested in any thoughts on this topic.
Best,
Chad