Here we will use the bayesdfa package to fit dynamic factor analysis (DFA) model to simulated time series data. In addition to working through an example of DFA for multivariate time series, we’ll apply bayesdfa routines for fitting hidden Markov models (HMMs) to the estimated trends to identify latent regimes. Most of the core functions of the package are included here, including `fit_dfa()`

and `find_dfa_trends`

for fitting DFA models, `plot_trends()`

, `plot_fitted()`

and `plot_loadings()`

for plotting estimates, `find_swans()`

for flagging extremes, `fit_regimes()`

and `find_regimes()`

for fitting HMM models, and `plot_regimes()`

for plotting HMM output.

Let’s load the necessary packages:

First, let’s simulate some data. We will use the built-in function `sim_dfa()`

, but normally you would start with your own data. We will simulate 20 data points from 4 time series, generated from 2 latent processes. For this first dataset, the data won’t include extremes, and loadings will be randomly assigned (default).

Next, we’ll fit a 1-trend, 2-trend, and 3-trend DFA model to the simulated time series using the `fit_dfa()`

function. Starting with the 1-trend model, we’ll estimate the posterior distributions of the trends and loadings. Note that this example uses 1 MCMC chain and 200 iterations — for real examples, you’ll want to use more (say 4 chains, 5000 iterations).

Convergence of DFA models can be evaluated with our `is_converged()`

function. This function takes a fitted object, and specified `threshold`

argument representing the maximum allowed Rhat value (default = 1.05). The convergence test isn’t that useful for a model with such a short number of iterations, but is called with

`## [1] TRUE`

Before we extract the trends from the model, we need to rotate the loadings matrix and trends. By default we use the varimax rotation, implemented in the `rotate_trends()`

function. An optional argument is the `conf_level`

argument, which calculates the specified confidence (credible) interval of the estimates (by default, this is set to 0.95).

The rotated object has several quantities of interest, including the mean values of the trends “trends_mean” and loadings “Z_rot_mean”,

```
## [1] "Z_rot" "trends" "Z_rot_mean" "Z_rot_median"
## [5] "trends_mean" "trends_median" "trends_lower" "trends_upper"
```

We can then plot the trends and intervals, with

We can also plot the estimated loadings (we’ll show that plot for the more complicated 2-trend model below because it’s not as interesting for the 1-trend model) and the fitted values. To plot the fitted values from the 1-trend model, we’ll use the `plot_fitted()`

function (predicted values can also be returned without a plot, with the `predicted()`

) function.

We can then plot the trends and intervals, faceting by time series, with

Moving to a more complex model, we’ll fit the 2-trend and 3-trend models. All other arguments stay the same as before,

```
f2 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, zscore = TRUE,
iter = 200, chains = 1, thin = 1
)
r2 <- rotate_trends(f2)
f3 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 3, zscore = TRUE,
iter = 200, chains = 1, thin = 1
)
r3 <- rotate_trends(f3)
```

The fits from the 2-trend model look considerably better than that from the 1-trend model,

The loadings from the 1-trend model aren’t as interesting because for a 1-trend model the loadings are a 1-dimensional vector. For the 2 trend model, there’s a separate loading of each time series on each trend,

```
## [,1] [,2]
## [1,] 0.78 -0.09
## [2,] -0.07 0.60
## [3,] 0.30 0.55
## [4,] -0.69 -0.26
```

These loadings can also be plotted with the `plot_loadings()`

function. This shows the distribution of the densities as violin plots, with color proportional to being different from 0.

Finally, we might be interested in comparing some measure of model selection across these models to identify whether the data support the 1-trend, 2-trend, or 3-trend models. The Leave One Out Information Criterion can be calculated with the `loo()`

function, for example the LOOIC for the 1-trend model can be accessed with

```
## Estimate SE
## elpd_loo -90.77104 5.221621
## p_loo 10.89671 1.303337
## looic 181.54208 10.443242
```

where 181.5420834 is the estimate and 10.4432421 is the standard error.

As an alternative to fitting each model individually as we did above, we also developed the `find_dfa_trends()`

to automate fitting a larger number of models. In addition to evaluating different trends, this function allows the user to optionally evaluate models with normal and Student-t process errors, and alternative variance structures (observation variance of time series being equal, or not). For example, to fit models with 1:5 trends, both Student-t and normal errors, and equal and unequal variances, the call would be

In this example, we’ll simulate data with an extreme anomaly. As before, this will be 20 data points from 4 time series, generated from 2 latent processes. The `sim_dfa()`

function’s arguments `extreme_value`

and `extreme_loc`

allow the user to specify the magnitude of the extreme (as an additive term in the random walk), and the location of the extreme (defaults to the midpoint of the time series). Here we’ll include an extreme value of 6,

Plotting the data shows the anomaly occurring between time step 9 and 10,

Though the plot is a little more clear if we standardize the time series first,

Instead of fitting a model with normal process deviations, we may be interested in fitting the model with Student-t deviations. We can turn on the estimation of `nu`

with the `estimate_nu`

argument. [Alternatively, nu can also be fixed a priori by setting the argument `nu_fixed`

]. Here’s the code for a 2-trend model with Student-t deviations,

```
t2 <- fit_dfa(
y = sim_dat$y_sim, num_trends = 2, zscore = TRUE,
iter = 200, chains = 1, thin = 1, estimate_nu = TRUE
)
```

Again we have to rotate the trends before plotting,

And the loadings,

One way to look for extremes is using the `find_swans()`

function, which evaluates the probability of observing a deviation in the estimated trend (or data) greater than what is expected from a normal distribution. This function takes a `threshold`

argument, which specifies the cutoff. For example, to find extremes greater than 1 in 1000 under a normal distribution, the function call is

Setting plot to TRUE also creates a time series plot that flags these values.

We can also look at the estimated `nu`

parameter, which shows some support for using the Student-t distribution (values greater than ~ 30 lead to similar behavior as a normal distribution),

```
## V1
## Min. : 5.807
## 1st Qu.:11.030
## Median :16.515
## Mean :19.342
## 3rd Qu.:24.094
## Max. :60.323
```