A motivating example

Suppose that a store is run and we want to perform a analysis of sales of a specific item. To set up the objective more concrete, we will put model under the sales-price relationship. Assumed that sales is composed of two components: an underlying level together with a determined amount of price. First, you may think that the sale is driven by a fixed trend (or a mean level) solely with some variations due to randomness. Then, you might come up with a simple model with no temporal structure: \[ y_{t} = \mu + \upsilon_{t}, \quad \upsilon_{t} \stackrel{i.i.d}{\sim} \mathcal{N}(0, V) \] An simple plot of this type is illustrated below:

library(rmdformats)
library(demography)
library(latex2exp)
fixed_lvl = 8
v = 2
y <- fixed_lvl + rnorm(100, 0, sd = sqrt(v))
plot(y, col = "red", type = "l", xlab = "time", ylab = TeX("y_t"))

lines(seq(1,100,1), rep(fixed_lvl, 100), col = "black")

legend("bottomleft", c("y_t", "level"), 
       lty = c(1,1), col = c("red", "green"))

The green line in the above plot is the expected value of \(y_{t}\) or our mean level, \(E(y_{t}) = \mu\) and the variance of \(y_{t}\) is \(V(y_{t}) = \upsilon\). If we want to capture the dynamics of the level, or we want to model some temporal structure of this series, we may write down the model as follows: \[\begin{align} y_{t} & = \mu_{t} + \upsilon_{t}, \quad \upsilon_{t} \stackrel{i.i.d}{\sim} \mathcal{N}(0, v), \\ \mu_{t} & = \mu_{t-1} + \omega_{t}, \quad \omega_{t} \stackrel{i.i.d}{\sim} \mathcal{N}(0, \omega) \end{align}\]

We could see that the equation depicting the evolution of the level is simply the random walk structure, or with evolution error \(\omega_{t}\). Additionally, The observational error \(\upsilon_{t}\) and the evolution error \(\omega_{t}\) are assumed to be mutally independent which both are normal random variables centering at \(0\). Hence, if we concern about the foregoing process of both equations, it is easy to see that for \(t = 1,2,\dots\), \[\begin{align} (Y_{t} \mid \mu_{t}) & \sim \mathcal{N}(\mu_{t}, V_{t}), \\ (\mu_{t} \mid \mu_{t - 1}) & \sim \mathcal{N}(\mu_{t-1}, W_{t}), \end{align}\]

We will given an example in the below figure

library(tibble)
library(dplyr)
library(ggplot2)
fixed_lvl_0 = 25
fixed_lvl <- matrix(0, ncol = 2, nrow = 101)
fixed_lvl[1,] <- fixed_lvl_0
v = 1
w = c(0.05, 0.5)
y = matrix(0,100, 2)
for (j in range(1,2)){
  for (i in 2:100){
    y[i-1,j] <- fixed_lvl[i-1,j] + rnorm(1, 0, sd = sqrt(v))
    fixed_lvl[i,j] <- fixed_lvl[i-1,j] + rnorm(1, 0, sd = sqrt(w[j]))
  }
}


y <- as_tibble(y)
fixed_lvl <- as_tibble(fixed_lvl)
y <- bind_cols(y, fixed_lvl[-1,])
colnames(y) <- c("W1","W2","Level1", "Level2")
time_ind <- seq(1,100,1)

Example series with W/V = 0.05

Example series with W/V = 0.5

DLM: An Introduction

Before diving into the general class of dynamics linear model, we will examine the first case of this class, namely first-order polynomial model or constant level plus random noise model. Even though this model is simple, it provides a fundamental baseline which will represent the evolution process of the DLMs.

Recalling the system of equations given in the above part:

\[\begin{align} y_{t} & = \mu_{t} + \upsilon_{t}, \quad \upsilon_{t} \stackrel{i.i.d}{\sim} \mathcal{N}(0, V), \\ \mu_{t} & = \mu_{t-1} + \omega_{t}, \quad \omega_{t} \stackrel{i.i.d}{\sim} \mathcal{N}(0, W) \end{align}\]

However, this time we denote that for all \(t\) and \(s\) with \(t \neq s\), \(\upsilon_{t}\) and \(\upsilon_{s}\), \(\omega_{t}\) and \(\omega_{s}\), and \(\omega_{t}\) and \(\upsilon_{s}\) are all independent. Moreover, it is assumed that the variances are dependent on time \(t\), \(V_{t}\) and \(W_{t}\), are known for each time \(t\). Observational and evolution equations are evolving at \(t = 1,2, \cdots,\) as below \[\begin{align} (Y_{t} &\mid \mu_{t}) \sim \mathcal{N}(\mu_t, V_t), \\ (\mu_{t} &\mid \mu_{t-1}) \sim \mathcal{N}(\mu_{t-1}, W_t), \end{align}\]

Two figures above illustrates chronological evolution of the level of series \(\{Y_{t}, \mu_{t}, t\}\). Starting value at \(\mu_{0} = 25\), and the observational variance \(V_{t} = 1\) is constant. Refering to the evolution variance, we consider two cases with \(W = 0.05\) and \(W = 0.5\). The former results in the ratio \(W/V = \frac{1}{20}\) and the latter gives \(W/V = \frac{1}{2}\). This ratio is also called signal-to-noise ratio and it reflects on the variation of structure of estimation and forecasting process. Further consideration or properties of this ratio will be provided in the next section.

Why we have to explore this kind of model.Apparently, this model often used in an reasonably effective way for depicting many real-world cases involving in short-term forecasting, especially in production planning and inventory control. For instance, we could imagine that the company want to predict the market demand for a specific product and \(\mu_{t}\) indicates the true underlying demand at time \(t\) with some random fluctuations \(\upsilon_{t}\), which stem from backorders. Then, in some local moments, the underlying demand \(\mu_{t}\) is specified as approximate constant. Moreover, the variations over a long horizon are expected as a pure stochastic process (since the evolution variances are zero mean processes), which implies that we are not interested in anticipating the form of this significant changes. Hence, one way to think \(\mu_{t}\) as a smoother \(\mu(t)\) with an associated Taylor series representation \[ \mu(t + \delta(t)) = \mu_{t} + O(t^2) \] where the higher-order terms are characterised as zero-mean noise.

## Forward filtering is completed!

## [1] "mt" "Ct" "at" "Rt" "ft" "Qt"

## Forecasting is completed!

## Forward filtering is completed!

## Forecasting is completed!