I’d have to know how all the pieces fit together to make forecasts. As an example which might be way off base, say I had time series data x_t which is determined by some exogenous process z_t, plus a local-linear trend:
\begin{align} x_t &= \beta z_t + \eta_t \\ \eta_t &= \eta_{t-1} + \alpha_{t-1} \\ \alpha_t &= \alpha_{t-1} + \zeta_t \end{align}
Where \eta_t is the (conditional) level of the series x_t, and \alpha_t is the “velocity” of change of the level. It’s a smooth trend because there’s no error term on \eta_t, only on \alpha_t. You could then define y_t = x_t + \xi_t as a measurement error process.
Anyway, If I was going to forecast this I couldn’t, because to know x_{t+h} I need to know z_{t+h}. If that first line of the equation wasn’t there, I could just iterate on random draws from \zeta_t \sim \text{Whatever} and forecast as far forward as I like.
If I had a secondary model for z_t, though, I could just directly fuse the two models together to make forecasts. The setup would be:
- Have the above model of x_t | z_t
- Estimate using data x_t, z_t, obtain posterior samples
idata_1 - Have some model of z_t.
- Write a 2nd pymc model which replaces z_t with data generated by my model from step (3), and recycles the posterior parameters obtained in step (2).
- Use
pm.sample_posterior_predictiveto combine the two models and obtain forecasts
So maybe this is a bit rambling, but the point was that because I started from the model that connected x_t and z_t , I could do the forecasts in a generative fashion. I guess in your case, x_t is sea level and z_t is antarctic ice-sheet. It doesn’t seem like you have z_t in your initial model. Nothing stops you from doing exactly what I’ve laid out here though without that modification – just estimate the time series x_t without z_t, then introduce z_t afterwards. That should be fine (as long as the posteriors you obtain on the x_t model are conditionally independent of z_t), as long as you can write down the relationship between x_t and z_t.