--- title: "farfima: Order Selection and Forecasting in Long-Memory Functional Time Series" author: "Chang Liu, Han Lin Shang, Garth Tarr and Samuel Müller" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{farfima: Order Selection and Forecasting in Long-Memory Functional Time Series} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} --- ```{r, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5, dpi = 200, out.width = "100%", eval = requireNamespace("demography", quietly = TRUE)) ``` ## 1. Introduction A functional time series is a sequence of curves observed over time, such as age-specific death rates in successive years. When the dependence between the curves decays slowly, the series has long memory. The functional autoregressive fractionally integrated moving-average (FARFIMA) model describes such series. The curves $X_t$ are fractionally differenced into short-memory curves $Y_t$, $$\nabla^d X_t(u) = Y_t(u), \qquad \nabla = 1 - B,$$ where $B$ is the backshift operator and $d$ is the memory parameter. The curves $Y_t$ follow a functional autoregressive moving-average (FARMA) model of orders $p$ and $q$, $$Y_t(u) - \sum_{i=1}^{p} (\Phi_i Y_{t-i})(u) = \eta_t(u) + \sum_{j=1}^{q} (\Theta_j \eta_{t-j})(u),$$ with autoregressive operators $\Phi_i$, moving-average operators $\Theta_j$ and innovations $\eta_t$. The farfima package fits this model in three steps. It estimates $d$, selects $p$ and $q$, and estimates the operators. It then forecasts the curves with conformal prediction bands. The main functions are: | Function | Purpose | |---|---| | `plot_curves()` | Rainbow plots of the curves | | `split_curves()` | Training, validation and test samples | | `farfima()` | Fit the model and calibrate the conformal prediction bands | | `fitted()`, `residuals()`, `plot()` | The in-sample fit | | `predict()` | Forecasts, and their evaluation on a test sample | | `plot()`, `as.data.frame()` | Plots and tables of the forecasts | | `estimate_memory()`, `frac_diff()`, `select_order()`, `estimate_operators()` | The steps of `farfima()` on their own | | `sim_farfima()` | Simulated curves | Each section below introduces one function and then applies it to French female mortality rates. ## 2. Data The data are French female death rates from the demography package, for the years 1901 to 2006. Ages 0 to 94 are kept as single years, and ages 95 and above are combined into one group, which avoids zero rates at the oldest ages. The rates are smoothed by weighted penalised regression splines, constrained to increase with age above 65, and taken on the log10 scale. In farfima, curves are stored as a matrix with grid points in rows and time in columns, here ages by years. `fts` objects from the rainbow and ftsa packages also work. ```{r, message = FALSE, warning = FALSE} library(farfima) library(demography) ``` ```{r} fr = extract.ages(extract.years(fr.mort, 1901:2006), 0:95) x = log10(smooth.demogdata(fr, method = "mspline", b = 65, weight = TRUE)$rate$female) dim(x) ``` ## 3. Rainbow plots: `plot_curves()` `plot_curves()` draws the curves coloured by time, from red for the earliest to violet for the latest. With `type = "rainbow"`, all curves are in one panel with a colour bar. With `type = "periods"`, each period has its own panel and the other curves are in grey. `periods` is the number of periods, or the first year of each period. The death rates fall over time at almost all ages. The curves of the war years stand out at young adult ages. ```{r, fig.height = 4} plot_curves(x, xlab = "Age", ylab = "log10 death rate") ``` ```{r, fig.height = 5.5} plot_curves(x, type = "periods", periods = c(1901, 1921, 1941, 1961, 1981), xlab = "Age", ylab = "log10 death rate") ``` ## 4. Training, validation and test samples: `split_curves()` `split_curves()` splits the curves in time order. `train` and `valid` are the proportions of curves in the training and validation samples, rounded down, and the remaining curves form the test sample. The model is estimated on the training sample, the conformal prediction bands are calibrated on the validation sample, and the forecasts are evaluated on the test sample. ```{r} blocks = split_curves(x, train = 0.6, valid = 0.2) sapply(blocks, ncol) ``` The training sample covers `r paste(range(colnames(blocks$train)), collapse = " to ")`, the validation sample `r paste(range(colnames(blocks$valid)), collapse = " to ")` and the test sample `r paste(range(colnames(blocks$test)), collapse = " to ")`. ## 5. Fitting the model: `farfima()` `farfima()` fits the model to the training sample in three steps. First, it estimates $d$ with `estimate_memory()`. Second, it centres and fractionally differences the training curves, drops the first `trim` differenced curves, and selects $p$ and $q$ as in `select_order()`. Third, it estimates the operators by the functional Hannan-Rissanen method, as in `estimate_operators()`. With a validation sample, it also forecasts up to `h` steps ahead from every validation origin, and uses the forecast errors to calibrate the conformal prediction bands for each horizon and level. The main arguments are: | Argument | Sets | Options | |---|---|---| | `d` | memory parameter | `"lw"`, `"dfa"`, `"lpwn1"`, `"lpwn2"`, `"bc_lpwn1"`, `"bc_lpwn2"`, or a number | | `order` | FARMA orders | `"ic_full"`, `"ic_auto"`, `"bj"`, `"cv"`, or `c(p, q)` | | `p`, `q` | candidate orders | e.g. `0:5` | | `method` | operator estimator | `"farma"`, or `"far"` for the FAR(1) benchmark | | `h` | longest horizon | `valid` needs at least `h + 1` curves | | `level` | levels of the conformal prediction bands | e.g. `c(0.8, 0.95)` | | `band` | conformal prediction bands | `"SR"`, `"NS"`, `"AC"`, `"PID"`, `"AcMCP"` | | `symmetric` | symmetric AC and PID bands | `TRUE` or `FALSE` | | `eta_grid` | candidate learning rates | e.g. `c(0, 0.005, 0.01, 0.02, 0.05, 0.1)` | | `scorecast` | scorecaster of the PID band | `"none"`, `"naive"`, `"theta"` | | `window` | window at each origin | `"expanding"` or `"fixed"` | | `trim` | differenced curves dropped | e.g. 20 | | `memory_args`, `order_args` | other settings of `estimate_memory()` and `select_order()` | e.g. `list(B = 50)` | The fit below estimates $d$ by BC-LPWN2, selects the orders by the functional information criterion over $p, q \in \{0, \dots, 5\}$, and calibrates the SR, NS and AC bands at the 80% and 95% levels for horizons up to 10 years. ```{r} fit = farfima(blocks$train, blocks$valid, d = "bc_lpwn2", order = "ic_full", p = 0:5, q = 0:5, h = 10, level = c(0.8, 0.95), band = c("SR", "NS", "AC"), trim = 20, seed = 1) fit ``` The estimate of $d$ is `r round(fit$d, 3)`. It is above 0.5, so the curves are non-stationary. The selected model for the differenced curves is FARMA(`r fit$order[1]`, `r fit$order[2]`). The fit is a list: ```{r} fit$d fit$order names(fit) head(fit$order_selection$table) head(fit$calibration) ``` ### Other settings Only the AC band: ```{r} fit_ac = farfima(blocks$train, blocks$valid, band = "AC", trim = 20, seed = 1, verbose = FALSE) fit_ac ``` Other estimators, known values and the FAR(1) benchmark: ```{r} fit_lw = farfima(blocks$train, blocks$valid, d = "lw", order = "bj", trim = 20, seed = 1, verbose = FALSE) fit_lw$d fit_lw$order fit_known = farfima(blocks$train, blocks$valid, d = 1, order = c(1, 0), trim = 20, verbose = FALSE) fit_known$order bench = farfima(blocks$train, blocks$valid, d = 1, method = "far", trim = 20, verbose = FALSE) bench ``` ### Without a validation sample Without `valid` there are no conformal prediction bands, and `order = "cv"` cannot be used. ```{r} fit_point = farfima(blocks$train, d = "lw", trim = 20, verbose = FALSE) fit_point ``` ```{r, error = TRUE} farfima(blocks$train, order = "cv", trim = 20, verbose = FALSE) ``` ## 6. In-sample fit: `fitted()`, `residuals()` and `plot()` Each fitted curve is the one-step-ahead forecast of a training curve from the curves before it, and each residual curve is the training curve minus its fitted curve. They cover the training curves left after trimming, and the first `max(p, q)` of them have no fit. `fitted()` gives the fitted curves on the scale of the data, or on the fractionally differenced scale with `scale = "differenced"`. `plot()` draws the fitted curves, or the residual curves with `type = "residuals"`. Here $\max(p, q) = `r max(fit$order)`$, so the first fitted curve, for `r colnames(fitted(fit))[1]`, is missing. `fitted(fit) + residuals(fit)` returns the training curves. ```{r} round(fitted(fit)[1:3, 1:4], 3) round(fitted(fit, scale = "differenced")[1:3, 1:4], 4) round(residuals(fit)[1:3, 1:4], 4) ``` The fitted curves follow the training curves closely. The largest residuals are in the early 1940s, during the Second World War. ```{r, fig.height = 4} plot(fit, xlab = "Age", ylab = "log10 death rate") ``` ```{r, fig.height = 4} plot(fit, type = "residuals", xlab = "Age") ``` ## 7. Forecasting: `predict()` `predict()` makes point forecasts up to `h` steps beyond the data, with the conformal prediction bands chosen in `farfima()`. Without `test`, the forecasts start after the last validation curve. With `test`, the forecasts are also made from every test origin and compared with the observed curves. The memory parameter and the orders stay fixed, and the operators and the mean are re-estimated at each origin. The evaluation has three tables. `point` gives the root mean squared prediction error (RMSPE) and the mean absolute prediction error (MAPE) of the point forecasts. `interval` gives the empirical coverage probability (ECP) of each band, the absolute difference between the ECP and the nominal level (CPD), and the mean interval score (MIS). `wis` gives the weighted interval score over the levels. Without `test`, `predict()` forecasts the `h` years after the end of the validation sample, 1984: ```{r} fc_future = predict(fit) fc_future ``` With `test`, it also forecasts from every test year and compares the forecasts with the observed curves: ```{r} fc = predict(fit, test = blocks$test) fc ``` ```{r} names(fc) names(fc$evaluation) round(fc$point_forecast[1:3, 1:4], 3) round(fc$bands$AC[["95%"]]$upper[1:3, 1:4], 3) fc$evaluation$point head(fc$evaluation$interval) head(fc$evaluation$wis) ``` `h` can be shorter than in the fit, and `band` can be a subset of the bands of the fit: ```{r} fc_ac5 = predict(fit, test = blocks$test, h = 5, band = "AC", verbose = FALSE) fc_ac5$evaluation$interval ``` A fit without a validation sample gives point forecasts only, and its test sample follows the training sample: ```{r} fc_point = predict(fit_point, test = cbind(blocks$valid, blocks$test), h = 5) fc_point ``` A test sample that does not follow on in time gives a warning: ```{r} fc_wrong = predict(fit_point, test = blocks$test, h = 5, verbose = FALSE) ``` ### AC, PID and AcMCP conformal prediction bands These bands are updated over the test sample from the error rates of earlier bands. With `symmetric = FALSE`, the lower and upper edges of the AC and PID bands are calibrated separately. `scorecast = "naive"` or `"theta"` adds a scorecaster to the PID band. The AcMCP band is always asymmetric. The learning rates are chosen on the validation sample: ```{r} fit_seq = farfima(blocks$train, blocks$valid, d = fit$d, order = fit$order, h = 10, band = c("NS", "AC", "PID", "AcMCP"), symmetric = FALSE, trim = 20, verbose = FALSE) fit_seq$calibration[fit_seq$calibration$h %in% c(1, 10), ] ``` ```{r, fig.width = 9, fig.height = 7} fc_seq = predict(fit_seq, test = blocks$test, verbose = FALSE) fc_seq plot(fc_seq, type = "horizons", h = c(1, 5, 10), xlab = "Age", ylab = "log10 death rate") ``` ## 8. Plotting forecasts: `plot()` `plot()` draws a forecast in four ways. `type = "curve"` draws the forecast `h` steps beyond the data with its conformal prediction band. With `target_test_curve = k`, it draws the forecast of the k-th test curve made `h` steps earlier, with the observed curve. `type = "horizons"` shows one column per band and one row per horizon. `type = "series"` shows a few grid points over time, and `type = "recent"` the last curves with the forecasts. The first plot is the one-year-ahead forecast for 2007 with its AC band. The second is the forecast of the 11th test year, 1995, made three years earlier, with the observed curve. ```{r, fig.height = 4} plot(fc, h = 1, band = "AC", xlab = "Age", ylab = "log10 death rate") ``` ```{r, fig.height = 4} plot(fc, h = 3, target_test_curve = 11, band = "AC", level = 0.95, xlab = "Age", ylab = "log10 death rate") ``` ```{r, fig.height = 7} plot(fc, type = "horizons", xlab = "Age", ylab = "log10 death rate") ``` ```{r, fig.height = 3.5} plot(fc, type = "series", at = c(0, 20, 50, 80), band = "AC", main = "Age", xlab = "Year", ylab = "log10 death rate") ``` ```{r, fig.height = 4.5} plot(fc, type = "recent", xlab = "Age", ylab = "log10 death rate") ``` ## 9. Tables: `as.data.frame()` `as.data.frame()` puts the forecasts in a table with one row per grid point, horizon, band and level. By default it gives the forecasts beyond the data. With `test = TRUE`, it gives the forecasts of the test curves, with their origins and the observed values. ```{r} head(as.data.frame(fc), 3) df_test = as.data.frame(fc, test = TRUE) head(df_test, 3) ``` The table makes summaries easy. For example, the coverage of the 95% AC band at `h = 1`: ```{r} ac95 = subset(df_test, band == "AC" & level == 0.95 & h == 1) mean(ac95$observed >= ac95$lower & ac95$observed <= ac95$upper) ``` ## 10. The steps on their own The three steps of `farfima()` can also be run on their own. The memory parameter is estimated on the original curves, and the orders and operators on short-memory curves. ### Memory parameter: `estimate_memory()` `estimate_memory()` estimates $d$ from the leading score process, the centred curves projected onto the leading eigenfunction of a covariance-type operator. The estimators are local Whittle (`"lw"`), detrended fluctuation analysis (`"dfa"`), local polynomial Whittle with noise (`"lpwn1"`, `"lpwn2"`), and their bias-corrected versions by the prefiltered sieve bootstrap (`"bc_lpwn1"`, `"bc_lpwn2"`). `m_const` sets the bandwidth, and `B` the number of bootstrap samples. With the same seed, BC-LPWN2 gives the estimate used by `fit`: ```{r} m = estimate_memory(blocks$train, method = "bc_lpwn2", seed = 1) m sapply(c("dfa", "lw", "lpwn1", "lpwn2"), function(meth) estimate_memory(blocks$train, method = meth)$d) ``` ### Fractional differencing: `frac_diff()` `frac_diff()` applies $\nabla^d = (1 - B)^d$ at every grid point. `farfima()` centres the training curves, fractionally differences them with `d` and drops the first `trim`. The same steps by hand give the short-memory curves `y`: ```{r} y = frac_diff(blocks$train - rowMeans(blocks$train), d = m$d)[, -(1:20)] dim(y) ``` A negative `d` undoes the differencing: ```{r} z = frac_diff(blocks$train, d = m$d) max(abs(frac_diff(z, d = -m$d) - blocks$train)) ``` ### FARMA orders: `select_order()` `select_order()` selects $p$ and $q$ for short-memory curves. `"ic_full"` minimises the adaptive functional Hannan-Quinn information criterion over all candidates, and `"ic_auto"` uses the same criterion with a stepwise search. `"bj"`, the functional Box-Jenkins method, tries the candidates by increasing $p + q$ and stops at the first whose residuals pass the functional portmanteau test. `"cv"` minimises the forecast error on a validation sample. `plot()` shows the criterion of each candidate, with the selected model boxed. ```{r} o_ic = select_order(y, method = "ic_full") o_ic select_order(y, method = "bj", seed = 1) ``` `"ic_auto"` warns when its largest model is unstable. A smaller range avoids this: ```{r} o_auto = select_order(y, method = "ic_auto") select_order(y, method = "ic_auto", p = 0:3, q = 0:3) ``` With `order = "cv"`, `farfima()` selects the orders by forecasting the validation curves on their original scale: ```{r} fit_cv = farfima(blocks$train, blocks$valid, d = m$d, order = "cv", h = 5, trim = 20, verbose = FALSE) fit_cv$order ``` ```{r, fig.width = 9, fig.height = 4} oldpar = par(mfrow = c(1, 2)) plot(o_ic) plot(fit_cv$order_selection) par(oldpar) ``` ### Operators: `estimate_operators()` `estimate_operators()` estimates the operators of a FARMA(p, q) model by the functional Hannan-Rissanen method. Each operator is stored as a grid by grid matrix, so applying it to a curve is a matrix product. `plot()` draws the operators as heatmaps, red for positive and blue for negative values. With the selected orders, it gives the operators of `fit`: ```{r} op = estimate_operators(y, order = o_ic) op all.equal(op[c("phi", "theta")], fit[c("phi", "theta")]) ``` ```{r, fig.width = 6, fig.height = 4.5} plot(op, xlab = "Age v", ylab = "Age u") ``` ## 11. Simulated data: `sim_farfima()` `sim_farfima()` simulates `n` curves from a FARFIMA(p, d, q) model on `no_grid` equally spaced points of [0, 1], with Brownian motion innovations. `ar_strength` and `ma_strength` set the operator norms. Simulated curves are useful to check the methods when the true model is known: ```{r} sim = sim_farfima(n = 500, d = 0.3, p = 1, q = 0, ar_strength = 0.5, seed = 1) sb = split_curves(sim) fit_sim = farfima(sb$train, sb$valid, d = "lpwn2", h = 5, trim = 50, verbose = FALSE) fit_sim$d fit_sim$order ``` ## 12. Unstable fits Large orders can give unstable operators. `farfima()` then warns, and forecast origins whose fit fails are skipped with a warning. ```{r} fit_33 = farfima(blocks$train, blocks$valid, d = fit$d, order = c(3, 3), h = 5, trim = 20, verbose = FALSE) fit_33 fc_33 = predict(fit_33, test = blocks$test, verbose = FALSE) ``` ## 13. Summary The farfima package fits FARFIMA models to functional time series in three steps, estimating the memory parameter, selecting the FARMA orders and estimating the operators, and forecasts the curves with conformal prediction bands. `farfima()` and `predict()` cover the usual workflow, and each step can also be run on its own. The help pages give the details of each function and its references.