Most of the theoretical papers on Variable Length Markov Chains
(VLMC, see `vignette("variable-length-markov-chains")`

) focus
on time series indexed by \(\mathbb{Z}\) which simplifies the analysis.
In practice time series are obviously finite which introduces some
difficulties associated to the initial observations. In particular, the
definition of a likelihood function is not entirely obvious.

This vignette discusses the case of VLMC, but the discussion applies to VLMC with covariates with minimal adaptation.

## Executive summary

For model selection, we recommend to use a likelihood function that
disregards the first few observations for which a (CO)VLMC cannot
provide a context. This is based on the *truncated* likelihood
function. For all practical uses such as prediction and sampling, we
recommend to use a notion of *extended* contexts.

## Likelihood functions for Markov chains

Let us consider a doubly infinite time series \(X_{i\in\mathbb{Z}}\) generated by some model \(\mathcal{M}\). The likelihood function associated to a finite observation of the time series, \((x_i)_{1\leq i\leq n}\), is \(\mathbb{P}_{\mathcal{M}}(X_1=x_1,\ldots,X_n=x_n)\).

If \(\mathcal{M}\) is a Markov chain of order \(d\), then we have \[ \mathbb{P}_{\mathcal{M}}(X_{d+1}=x_{d+1},\ldots,X_n=x_n)=\prod_{i=d+1}^n\mathbb{P}_{\mathcal{M}}(X_i=x_i|X_{i-1}=x_{i-1},\ldots,X_{i-d}=x_{i-d}). \] Thus if we know only \(x_{1\leq i\leq n}\) we can compute only the likelihood of \((x_i)_{d+1\leq i\leq n}\).

In practice, comparing models with different orders should be done only using the likelihood functions based on the same subset of the observed time series, i.e. using the highest order. As pointed out in several papers this has no impact on asymptotic results (see e.g. Garivier, A. (2006), Consistency of the unlimited BIC context tree estimator. IEEE Transactions on Information Theory, 52 (10) 4630–4635).

## Likelihood functions for VLMC

### Truncated solution

The simplest way to define a likelihood function for a single VLMC is to consider it as a Markov chain of the order given by the length of its longest context (i.e. of the order of the VLMC).

### Specific contexts

Another approach considers the fact that the past of \((x_i)_{1\leq i\leq n}\) is unknown and replace it by a collection of specific contexts that summarize this unknown past. This is used in Garivier’s paper cited above. By definition, each observation that does not have a actual context in the VLMC appears only once and thus is perfectly predicted by the empirical distribution associated to it. It has therefore a likelihood of 1. In practice, this amounts to identifying \(\mathbb{P}_{\mathcal{M}}(X_1=x_1,\ldots,X_n=x_n)\) to \(\mathbb{P}_{\mathcal{M}}(X_{d+1}=x_{d+1},\ldots,X_n=x_n)\).

In terms of parameters for AIC/BIC calculation, this corresponds to adding a specific parameter for each of the \(d\) initial values, where \(d\) is the order of the VLMC. Notice that this departs from Garivier’s approach cited above (in this paper, each initial value is associated to a full context and thus to \(|S|-1\) parameters if \(S\) is the state space).

### Extended contexts

Another approach considers extended/approximate contexts for the \(d\) initial values. Indeed each observation \(x_i\) for \(1\leq i\leq d\) can be considered has having its context partially determined by \((x_1,\ldots, x_{i-1})\), in particular by the empty context for \(i=1\). Let us consider for instance \(x_1\) and the empty context (the root of the context tree). If \(d\geq 1\), we cannot determine the context of \(x_1\) without values of \((x_{-d+1},\ldots, x_0)\) even if the empty context is a valid one. However, we can assign this empty context to \(x_1\) because of the lack of information. More generally, we can traverse the context tree using as many past values as available and stop in the corresponding node which is then interpreted as a context using the frequencies collected during the construction of the VLMC.

In terms of parameters, this adds to the VLMC an additional extended context for each node of the context tree which is not a context. For instance, in a binary state space \(S=\{0, 1\}\), one may consider a Markov Chain of order 1. In this case, we have two contexts \(0\) and \(1\), and the root node of the context tree is not a proper (empty) context. To compute the likelihood of the first observation, we need therefore a new extended context, the empty one.

### Complete example

Let us revisit the California earth quakes example proposed in
`vignette("variable-length-markov-chains")`

. The model is
obtained as follows (see the vignette for details):

```
California_centre <- data.frame(longitude = -119.449444, latitude = 37.166111)
distances <- geodist(globalearthquake[, c("longitude", "latitude")],
California_centre,
measure = "geodesic"
)
California_earth_quakes <- globalearthquake[distances < 2e6, ] ## distances are in meters
California_weeks <- rep(0, max(globalearthquake$nbweeks))
California_weeks[California_earth_quakes$nbweeks] <- 1
California_weeks_earth_quakes_model <- tune_vlmc(California_weeks,
initial = "truncated",
save = "all"
)
model <- as_vlmc(California_weeks_earth_quakes_model)
draw(model, prob = FALSE)
#> * (5126, 1291)
#> +-- 0 (4185, 940)
#> | +-- 0 (3460, 725)
#> | | '-- 0 (2883, 577)
#> | | '-- 0 (2421, 462)
#> | | '-- 1 (358, 103)
#> | '-- 1 (724, 215)
#> '-- 1 (940, 351)
```

The optimal model according to the BIC has an order of 5 and thus the simple truncated log likelihood is obtained by considering only the observations starting at index 6. We disregard the first 5 observations. The corresponding log likelihood is

```
loglikelihood(model, initial = "truncated")
#> 'log Lik.' -3182.077 (df= 6, nb obs.= 6412, initial="truncated")
```

Using specific contexts amounts to assuming perfect predictions for the first five observations. The log likelihood is not modified, but it covers now the full time series, i.e. 6417 observations. The number of parameters is increased as the initial value is now associated to a specific parameter leading to a total of 11 parameters.

```
loglikelihood(model, initial = "specific")
#> 'log Lik.' -3182.077 (df= 11, nb obs.= 6417, initial="specific")
```

The extended context approach is the most complex. For the first observation, we use the empty context, i.e. the root of the context tree. The associated empirical distribution is \(\mathbb{P}(X_1=1)=\frac{1291}{5126+1291}\simeq0.201\). Its contribution to the log likelihood is therefore \(\log \mathbb{P}(X_1=0)\simeq -0.225\) as the first observation is equal to 0.

As the root node is not a proper context, the specification of its associated empirical distribution contributes to the total number of parameters of the model (i.e. it adds a parameter to the total).

For the second observation, \(X_2=0\), the candidate context is \(0\). However, \(0\) is again not a proper context we should normally look for \(X_{0}\) and older values to find a proper context. In the extended approach, we consider the empirical distribution of values following a 0 in the time series, which is given by \(\mathbb{P}(X_t=1|X_{t-1}=0)=\frac{940}{4185+940}\simeq0.183\). The contribution of this observation to the log likelihood is therefore \(\log \mathbb{P}(X_t=0|X_{t-1}=0)\simeq -0.203\). In addition, this extended context adds again a parameter to the total.

The following observations \(X_3\), \(X_4\) and \(X_5\) have all proper contexts and contribute in a normal way to the (log) likelihood without the need for additional parameters. Notice that the number of additional parameters does not depend on the initial sequence but on the structure of the context tree. Notice also that while the corresponding nodes are contexts, they are not the contexts of those three values. Those contexts cannot be computed due to the lack of older values.

The final value is

```
loglikelihood(model, initial = "extended")
#> 'log Lik.' -3183.052 (df= 8, nb obs.= 6417, initial="extended")
```

## Monotonicity

### Nested models

For a given time series, the candidate VLMCs generated by the context algorithm follow a nested structure associated to the pruning operation: the most complete context tree is pruned recursively in order to generate less and less complex trees. A context in a small tree is also a suffix of context of a larger tree. The inclusion order is total unless some cut off values are identical.

A natural expected property of likelihood functions is to observe a decrease in likelihood for a given time series when switching from a model \(m_1\) to a simpler model \(m_2\) provided their parameters are estimated by maximum likelihood using this time series, in particular when \(m_2\) is nested in \(m_1\) in the sense of the likelihood-ratio test.

### Theoretical analysis

Let us consider the case of two VLMC models \(m_1\) and \(m_2\) where \(m_2\) is obtained by pruning \(m_1\) (both estimated on \((x_i)_{1\leq i\leq n}\)). Let us first consider \((x_i)_{i>d}\) where \(d\) is the order of \(m_1\). The contexts of those values are well defined, both in \(m_1\) and \(m_2\). The corresponding probabilities can be factorised according to the contexts as follows (for \(m_1\)) \[ \mathbb{P}_{m_1}(X_{d+1}=x_{d+1},\ldots,X_n=x_n)=\prod_{c\in m_1}\prod_{k,\ ctx(m_1, x_k)=c}\mathbb{P}_{m_1}(X_k=x_k|ctx(m_1, x_k)=c), \] where with use \(ctx(m_1, x_k)\) to denote the context of \(x_k\) in \(m_1\) and \(c\in m_1\) to denote all contexts in \(m_1\). We have obviously a similar equation for \(m_2\).

As \(m_2\) is included in \(m_1\) we know that each \(c\in m_2\) is also the suffix of a context of \(m_1\). When \(c\) is a context in both models, they concern the same subset of the time series and use therefore the same estimated conditional probabilities, leading to identical values of \[\mathbb{P}_{m_1}(X_k=x_k|ctx(m_1, x_k)=c)=\mathbb{P}_{m2}(X_k=x_k|ctx(m_2, x_k)=c).\]

When \(c\) is the suffix of a context in \(m_1\), there is a collection of contexts \(c'\) for which \(c\) is also a suffix. We have then \[ \{k\mid ctx(m_2, x_k)=c\}=\bigcup_{c'\in m_1, c\text{ is a suffix of }c'}\{j\mid ctx(m_1, x_j)=c'\}. \] In words, the collection of observations whose context in \(m_1\) has \(c\) as a suffix is equal to the collection of observations chose context in \(m_2\) is \(c\). Because the conditional probabilities are estimated by maximum likelihood, we have then \[ \begin{multline*} \prod_{\{k\mid ctx(m_2, x_k)=c\}}\mathbb{P}_{m_2}(X_k=x_k|ctx(m_2, x_k)=c)\leq\\ \prod_{\{l\mid ctx(m_1, x_k)=c', c\text{ is a suffix of }c'\}}\mathbb{P}_{m_1}(X_l=x_l|ctx(m_1, x_k)=c'). \end{multline*} \] Thus overall, as expected, \[ \mathbb{P}_{m_2}(X_{d+1}=x_{d+1},\ldots,X_n=x_n)\leq \mathbb{P}_{m_1}(X_{d+1}=x_{d+1},\ldots,X_n=x_n). \] However, the models may not have the same order. Fortunately, as probabilities are not larger than 1, we have also \[ \mathbb{P}_{m_2}(X_{d_{m_2}+1}=x_{d_{m_2}+1},\ldots,X_n=x_n)\leq \mathbb{P}_{m_1}(X_{d_{m_1}+1}=x_{d_{m_1}+1},\ldots,X_n=x_n), \] where \(d_m\) denotes the order of model \(m\).

In conclusion, likelihood functions based on *truncation* or
on *specific* contexts are non increasing when one moves from a
VLMC to one of its pruned version.

The case of *extended* contexts is more complex. Using the
hypotheses as above, the extended likelihood includes approximate
contexts for observations \(x_{1},\ldots,
x_{d_{m_1}}\) and for \(x_{1},\ldots,
x_{d_{m_2}}\). Let us consider the case where \(d_{m_1}=d_{m_2}+1\) (without loss of
generality). The only difference in the extended likelihoods is then
\(\mathbb{P}_{m_1}(X_{d_{m_1}}=x_{d_{m_1}}|ectx(m_1,
x_{d_{m_1}}))\) which is computed using the extended context
\(ectx(m_1, X_{d_{m_1}})\) and \(\mathbb{P}_{m_2}(X_{d_{m_1}}=x_{d_{m_1}}|ctx(m_2,
x_{d_{m_1}}))\) which is computed using the true context.

Notice that \[
(x_1,\ldots,x_{d_{m_1}-1},x_{d_{m_1}})=(x_1,\ldots,x_{d_{m_2}-1},x_{d_{m_2}},
x_{d_{m_1}})
\] Thus one computes the (extended) context of \(x_{d_{m_1}}\), the \(d_{m_2}\) first steps are obviously
identical in \(m_1\) and in \(m_2\), as \(m_2\) has been obtained by pruning \(m_1\). Depending on the structure of the
context tree, it may be possible possible to determine the true of
context of \(x_{d_{m_1}}\) in both
trees and thus the corresponding probabilities will be identical. The
only non obvious situation is when the context of \(x_{d_{m_1}}\) cannot be determine in \(m_1\). This is only possible if the context
in \(m_2\) is of length \(d_{m_2}\) and the corresponding leaf was an
internal node in \(m_1\). Then we use
as the *extended* context in \(m_1\) the *normal* context of \(m_2\) and therefore \[
\mathbb{P}_{m_1}(X_{d_{m_1}}=x_{d_{m_1}}|ectx(m_1,
x_{d_{m_1}}))=\mathbb{P}_{m_2}(X_{d_{m_1}}=x_{d_{m_1}}|ctx(m_2,
x_{d_{m_1}})).
\]

Thus in all cases, in the *extended* context interpretation,
moving from an extended context to a true context does not change the
probability included in the likelihood. Therefore, the extended
likelihood is also non increasing with the pruning operation.

### Experimental illustration

We used the saving options of `tune_vlmc()`

to keep all
the models considered during the model selection process in the
California earth quakes analysis. We can use them to illustrate non
decreasing behaviour of the likelihoods. For the *extended*
likelikood, we have:

```
CA_models <- c(
list(California_weeks_earth_quakes_model$saved_models$initial),
California_weeks_earth_quakes_model$saved_models$all
)
CA_extended <- data.frame(
cutoff = sapply(CA_models, \(x) x$cutoff),
loglikelihood = sapply(CA_models, loglikelihood,
initial = "extended"
)
)
ggplot(CA_extended, aes(cutoff, loglikelihood)) +
geom_line() +
xlab("Cut off (native scale)") +
ylab("Log likelihood") +
ggtitle("Extended log likelihood")
```

For the *specific* likelihood we have:

```
CA_specific <- data.frame(
cutoff = sapply(CA_models, \(x) x$cutoff),
loglikelihood = sapply(CA_models, loglikelihood,
initial = "specific"
)
)
ggplot(CA_specific, aes(cutoff, loglikelihood)) +
geom_line() +
xlab("Cut off (native scale)") +
ylab("Log likelihood") +
ggtitle("Specific log likelihood")
```

Notice that the time series contains 6417 observations and the maximum
order considered by `tune_vlmc()`

is 29, thus we do not
expect to observe large differences between the different log
likelihoods. This is illustrated on the following figure:

```
CW_combined <- rbind(
CA_extended[c("cutoff", "loglikelihood")],
CA_specific[c("cutoff", "loglikelihood")]
)
CW_combined[["Likelihood function"]] <- rep(c("extended", "specific"), times = rep(nrow(California_weeks_earth_quakes_model$results), 2))
ggplot(
CW_combined,
aes(cutoff, loglikelihood, color = `Likelihood function`)
) +
geom_line() +
xlab("Cut off (native scale)") +
ylab("Log likelihood") +
ggtitle("Log likelihood")
```

The case of the *truncated* likelihood is more complex. As
noted above, the numerical values of the *truncated* likelihood
are identical to the values of the *specific* likelihood and thus
the above graphical representations are also valid for it. However care
must be exercised when using the *truncated* likelihood for model
selection, as explained below.

## Model selection

Optimal (CO)VLMC models are generally selected via penalized
likelihood approaches, with a preference for the BIC, based on its
asymptotic consistency. A natural question is to what extent the
different likelihood functions proposed above are adapted for model
selection when combined with a penalty. Notice that consistency results
for the BIC are generally obtained with a *truncated* likelihood
function.

### Specific and extended likelihood

The *specific* and the *extended* likelihood functions
do not introduce any obvious difficulty. In particular, they work with
the full data set as they extend the VLMC model with specific/extended
contexts for the initial values. However, they tend to penalize complex
models more than the *truncated* likelihood. For instance, the
model selected on the California earth quakes is simpler than the one
selected with the *truncated* likelihood (as well as by the
*specific* one):

```
CA_model_extented <- tune_vlmc(California_weeks, initial = "extended")
model_extended <- as_vlmc(CA_model_extented)
draw(model_extended, prob = FALSE)
#> * (5126, 1291)
#> '-- 1 (940, 351)
```

This is also the case for, e.g., the sun spot time series used in the introduction to the package, as shown below:

```
sun_activity <- as.factor(ifelse(sunspot.year >= median(sunspot.year), "high", "low"))
sun_model_tune_truncated <- tune_vlmc(sun_activity, initial = "truncated")
draw(as_vlmc(sun_model_tune_truncated))
#> * (0.5052, 0.4948)
#> +-- high (0.8207, 0.1793)
#> | +-- high (0.7899, 0.2101)
#> | | +-- high (0.7447, 0.2553)
#> | | | +-- high (0.6571, 0.3429)
#> | | | | '-- low (0.9167, 0.08333)
#> | | | '-- low (1, 0)
#> | | '-- low (0.96, 0.04)
#> | '-- low (0.9615, 0.03846)
#> '-- low (0.1888, 0.8112)
#> +-- high (0, 1)
#> '-- low (0.2328, 0.7672)
#> +-- high (0, 1)
#> '-- low (0.3034, 0.6966)
#> '-- high (0.07692, 0.9231)
```

```
sun_model_tune_extended <- tune_vlmc(sun_activity, initial = "extended")
draw(as_vlmc(sun_model_tune_extended))
#> * (0.5052, 0.4948)
#> +-- high (0.8207, 0.1793)
#> | '-- high (0.7899, 0.2101)
#> | '-- high (0.7447, 0.2553)
#> | +-- high (0.6571, 0.3429)
#> | | '-- low (0.9167, 0.08333)
#> | '-- low (1, 0)
#> '-- low (0.1888, 0.8112)
#> +-- high (0, 1)
#> '-- low (0.2328, 0.7672)
#> '-- high (0, 1)
```

In this latter case, the *specific* likelihood gives the same
results as the *truncated* one.

We observe a similar behaviour on a simple second order Markov chain generated as follows:

```
TM0 <- matrix(c(0.7, 0.3, 0.4, 0.6),
ncol = 2,
byrow = TRUE
)
TM1 <- matrix(c(0.4, 0.6, 0.8, 0.2),
ncol = 2,
byrow = TRUE
)
init <- c(0, 1)
dts <- c(init, rep(NA, 500))
set.seed(0)
for (i in 3:length(dts)) {
if (dts[i - 1] == 0) {
probs <- TM0[dts[i - 2] + 1, ]
} else {
probs <- TM1[dts[i - 2] + 1, ]
}
dts[i] <- sample(0:1, 1, prob = probs)
}
```

Once again, the *extended* likelihood tends to over penalize
“complex” models

```
MC_extended <- tune_vlmc(dts, initial = "extended", save = "all")
draw(as_vlmc(MC_extended))
#> * (0.5458, 0.4542)
#> '-- 1 (0.5746, 0.4254)
#> +-- 0 (0.374, 0.626)
#> '-- 1 (0.8454, 0.1546)
```

while thishappens neither for the *truncated* likelihood

```
MC_truncated <- tune_vlmc(dts, initial = "truncated", save = "all")
draw(as_vlmc(MC_truncated))
#> * (0.5458, 0.4542)
#> +-- 0 (0.5201, 0.4799)
#> | +-- 0 (0.6408, 0.3592)
#> | '-- 1 (0.3923, 0.6077)
#> '-- 1 (0.5746, 0.4254)
#> +-- 0 (0.374, 0.626)
#> '-- 1 (0.8454, 0.1546)
```

no for the *specific* likelihood

```
MC_specific <- tune_vlmc(dts, initial = "specific", save = "all")
draw(as_vlmc(MC_specific))
#> * (0.5458, 0.4542)
#> +-- 0 (0.5201, 0.4799)
#> | +-- 0 (0.6408, 0.3592)
#> | '-- 1 (0.3923, 0.6077)
#> '-- 1 (0.5746, 0.4254)
#> +-- 0 (0.374, 0.626)
#> '-- 1 (0.8454, 0.1546)
```

If we increase the number of observations in the synthetic example,
to e.g. 5000, the three likelihood functions conduct to the same optimal
(and true) model, as expected. On smaller data sets, the use of the
*truncated* seem to be more adapted.

### Truncated likelihood

However, the *truncated* likelihood function is problematic
when used naively. The difficulty comes from the discarded observations:
two VLMC models with different orders are evaluated on different data
sets. If we consider for instance the BIC as an approximation of the
logarithm of the evidence of the data given the model, it is obvious
that one cannot compare directly two models on different data sets.

A possible solution for model selection based on the
*truncated* likelihood consists in choosing a maximal order, say
\(D\) and in evaluating the models only
on \((x_i)_{D+1\leq i\leq n}\) so that
contexts are always computable. This amounts to additional truncation
for simpler models. The solution is used by `tune_vlmc()`

and
`tune_covlmc()`

. All the examples given above have been
constructed using this approach.

## Coherence with other uses of a VLMC

### Sampling

As detailed in `vignette("sampling")`

a VLMC model can be
used to generate new discrete time series based on the conditional
probability distributions associated to the contexts. However, raw VLMC
models do not specify distributions for the initial values for which no
proper context exists.

The difficulty is generally circumvented by using a constant
initialisation coupled with a burn in phase. This is supported by the
theoretical results on VLMC bootstrap proved by Bühlmann and Wyner in their seminal
paper. Indeed the hypothesis used in the paper ensure an exponential
mix-in for Markov Chain and thus the initial values play essentially no
role in the stationary distribution, provided the burn in period is
“*long enough*”.

In practice, it is difficult to verify that the conditions of the
theorem apply and to turns them into actual numerical values of a
*long enough* burn in time. Thus we use in
`simulate.vlmc()`

the extended contexts proposed above. This
extends the VLMC into a full model, with specific probability
distributions for the initial values. This does not prevent e.g. slow
mix-in and this does not change the stationary distribution when it
exists, but this gives some coherence between the *extended*
likelihood function and sampling. Notice that we also support burn in
period as well as a manual specification of the initial values of a
sample.

### Prediction

A VLMC can also be used for one step ahead prediction of a time
series (or even for multiple steps ahead), as implement in
`predict.vlmc()`

. This poses obviously the same problem as
sampling or likelihood calculation for the initial values. We use in
`predict.vlmc()`

the extended contexts proposed above, using
the same extended VLMC model principle. This provides full coherence
between sampling, likelihood calculation and prediction.

The `metrics.vlmc()`

function computes predictive
performances of a VLMC on the time series used to estimate it.
Predictions used for those computations are also based on the extended
context principle.