After algebraic gymnastics with the backshift operator $\text{B}$ (i.e., $\text{B}y_t=y_{t - 1}$) I thought I found a convenient dynamic representation for a nonlinear model, but the representation fails in simulation. The situation is described below, but the high level questions are:

  • Are you allowed to multiply a transfer function numerator and denominator by $\text{B}$, leading to higher-order lags?
  • If your shift operator algebra leads to a problematic error structure for OLS or GLS estimation, what are your options?

Situation: In an athletic performance context, assume fitness $x_t$ and fatigue $w_t$ decay by proportions $\theta_x$ and $\theta_w$ per unit of time and are driven by exogenous training stimulus $u_t$:

$$ x_t = \theta_x x_{t - 1} + \theta_x u_{t - 1} \\ w_t = \theta_w w_{t - 1} + \theta_w u_{t - 1}. $$ Performance is modeled to be the effective net of fitness and fatigue plus white noise:

$$ y_t = \mu + k_x x_t + k_w w_t + \epsilon_t, $$ where $\epsilon_t \sim N(0, \sigma^2)$. Using the backshift operator $\text{B}$, this can be reframed as: $$y_t - \mu = \left[\frac{k_x\theta_x\text{B}}{1 - \theta_x\text{B}} + \frac{k_w\theta_w\text{B}}{1 - \theta_w\text{B}}\right] u_t + \epsilon_t.$$

I tried to combine the fractions in the brackets by getting both denominators to $(1 - \theta_x\text{B})(1 - \theta_w\text{B})$, or $$1 - (\theta_x + \theta_w) \text{B} + \theta_x \theta_w \text{B}^2.$$ This seemed great, because I could multiply by both sides of the equation by this "quantity" to knock out the denominator components. What's left is apparently a regression model of $y_t$ vs $y_{t-1}, y_{t-2}, u_{t-1}, u_{t-2}$ with MA(2) errors.

However, when I tried this with simulated data, I was unable to recover the theoretical coeffients of the dynamic model. I later realized that -($\theta_x + \theta_w$), one of the MA(2) parameters and the theoretical coefficient of $y_{t - 1}$, is likely to be absolutely greater than 1. Did I push this method too far? What rules did I break?

Update To follow up on @whuber's suggestion that there might be something wrong with the model, allow me to provide a bit more context. I've been studying the fitness-fatigue model of athletic performance, and I wrote this blog article which contains code to simulate data and recover parameters from the following model: enter image description here Later I came across an article from Kolossa et al where they put the fitness-fatigue model in the Kalman Filter framework. This is how they structured it, where A is the transition matrix, B multiplies the control (training) inputs, and C multiplies the latent vector to get to the observable performance: enter image description here.

I'm using a slightly simpler version where the bottom entry of the B-matrix is just like the top and I have both entries of C as free parameters, but I did test it out and it seemed to work. This dynamic representation was what motivated me to see if there was an even simpler dynamic regression implementation.

  • $\begingroup$ You seem to be using "$B$" to represent three distinct operators: one of them backshifts $w,$ another backshifts $u,$ and another backshifts $(w,u).$ $\endgroup$ – whuber Apr 11 '19 at 14:30
  • $\begingroup$ Is that not okay? I'm just thinking of d/dt in the continuous case that could be applied to x(t), y(t), so on and so forth. $\endgroup$ – Ben Ogorek Apr 11 '19 at 14:33
  • $\begingroup$ Even with derivatives that would be a mistake, because it would be tantamount to confusing $\partial/\partial t$ with $\partial/\partial s$ applied to a two-variable function $f(s,t).$ $\endgroup$ – whuber Apr 11 '19 at 14:43
  • $\begingroup$ Ok I need to think about this a bit. If you have time to write it up as an answer, I'd be indebted! $\endgroup$ – Ben Ogorek Apr 11 '19 at 14:54

I wonder about the model.

Here's why. Let's assume (as is implied) that $\theta_w$ and $\theta_x$ are nonzero. Notice that

$$\theta_w x_t - \theta_x w_t = \theta_w(\theta_x x_{t-1} + \theta_x u_{t-1}) - \theta_x(\theta_w w_{t-1} + \theta_w u_{t-1}) = 0.$$

Thus, you only need to keep track of one variable--say $w_t$--and you can reconstruct the other as

$$x_t = \frac{\theta_x}{\theta_w} w_t.$$

Consequently, setting $\kappa = k_x \theta_x/\theta_w + k_w,$

$$y_t = k + k_x x_t + k_w w_t + \epsilon_t = k + \kappa\,w_t + \epsilon _t$$

reduces this to a problem in which $x_t$ isn't involved.

It doesn't seem worthwhile proceeding with any analysis until we can resolve whether the model itself expresses your objectives correctly.

  • $\begingroup$ It's a real model, though I could have made a mistake expressing it above. This article provides code to simulate it and recover parameters using nonlinear least squares: towardsdatascience.com/…. For the dynamic representation above, I used this article's kalman filter form: degruyter.com/downloadpdf/j/ijcss.2017.16.issue-2/… (Equation 6). I will add more information about the original model form above in case that's where the error is. $\endgroup$ – Ben Ogorek Apr 11 '19 at 17:08
  • $\begingroup$ I couldn't find your model in either reference. I suspect you might have misunderstood their notation. $\endgroup$ – whuber Apr 11 '19 at 20:17

After some research, I have an answer and a simulation. A very helpful paper was Identification of Multiple-Input Transfer Function Models by Liu & Hanssens (1982). In the article they discuss the usefulness of a "common filter," e.g., premultipling right and left-hand side by $(1-\rho \text{B})$ where $\text{B}$ is the backshift operator, to get rid of autoregressive components with roots close to one (primarily for numerical reasons, they say).

A key point is that a common filter does not change the transfer function weights, so you can apply multiple common filters in succession, and besides causing complexity and losing some rows due to lags, you have not destroyed the relationship between the input and output series.

While I didn't come across anything that suggested that you can't treat the backshift operator like an algebraic quantity, what does it mean to say divide both sides of the equation by just $\text{B}$, for example? (Nothing that I can think of.) I do believe algebraic mistakes were the cause of some (but not all) of my confusion earlier, but in any case, when I started working in terms of common filters, the algebra became easier.


Part of the simulation is outsourced to this gist but it is simulating the model described in the original post. Here is the remaining simulation, where the user can replicate by running the gist or changing the path in read.csv. $\theta_x$, $\theta_w, k_x, k_w$ are theta1, theta2, k1 and k2 below, and the theoretical values come from coefficients after doing the common filter algebra.


# Regression parameters:

mu <- 496
theta1 <- exp(-1 / 60)
theta2 <- exp(-1 / 13)
k1 <- .07
k2 <- -.27
# Theoretical coefficient for intercept is 0.6069828
-1. * mu * (-1 + theta1 + theta2 - theta1 * theta2)

# Theoretical coefficient for performance lagged once is 1.909433
theta1 + theta2

# Theoretical coefficient for performance lagged twice is -0.9106563
-1. * theta1 * theta2

# Theoretical coefficient for training lagged once is -0.1811665
k1 * theta1 + k2 * theta2

# Theoretical coefficient for training lagged twice is 0.1821313
-1. * theta1 * theta2 * (k2 + k1)

# run code here: https://gist.github.com/baogorek/6d682e42079005b3bde951e98ebae89e
# or get file here: https://drive.google.com/open?id=1kk40wiVYzPXOkrPffU55Vzy-LLTrgAVh
train_df <- read.csv("/mnt/c/devl/data/train_df.csv")

train_aug <- train_df %>%
  mutate(perf_lag1 = lag(perf, n = 1, order_by = day),
         perf_lag2 = lag(perf, n = 2, order_by = day),
         train_lag1 = lag(w, n = 1, order_by = day),
         train_lag2 = lag(w, n = 2, order_by = day))

my_gls <- gls(perf ~ perf_lag1 + perf_lag2 + train_lag1 + train_lag2,
              data = train_aug[3:nrow(train_aug), ],
              corARMA(form = ~day, p = 0, q = 2))

my_lm <- lm(perf ~ perf_lag1 + perf_lag2 + train_lag1 + train_lag2,
            data = train_aug[3:nrow(train_aug), ])

Now check out the coefficient estimates from the GLS regression:

Correlation Structure: ARMA(0,2)
Formula: ~day
Parameter estimate(s):
Theta1     Theta2
-1.9059497  0.9117409

                 Value  Std.Error    t-value p-value
(Intercept)  0.6571088 0.11700730    5.61596       0
perf_lag1    1.9187158 0.00815689  235.22646       0
perf_lag2   -0.9200058 0.00815495 -112.81568       0
train_lag1  -0.1662026 0.02238219   -7.42566       0
train_lag2   0.1664704 0.02241510    7.42671       0

which clearly recovers the true values, and the OLS regression:

            Estimate Std. Error t value Pr(>|t|)
(Intercept) 33.99821   12.09470   2.811 0.005327 **
perf_lag1    0.48066    0.05619   8.553 1.18e-15 ***
perf_lag2    0.46189    0.05504   8.393 3.45e-15 ***
train_lag1  -0.15602    0.04406  -3.541 0.000475 ***
train_lag2  -0.02346    0.04516  -0.520 0.603807

which is awful! This is an important point: our common filter induced a gnarly MA(2) error structure on the response, and this is apparently enough to really throw off the estimates. I guess I don't understand why there is bias instead of simply more variability, so I'd appreciate any insight there.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.