## ----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))

## ----message = FALSE, warning = FALSE-----------------------------------------
library(farfima)
library(demography)

## -----------------------------------------------------------------------------
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)

## ----fig.height = 4-----------------------------------------------------------
plot_curves(x, xlab = "Age", ylab = "log10 death rate")

## ----fig.height = 5.5---------------------------------------------------------
plot_curves(x, type = "periods", periods = c(1901, 1921, 1941, 1961, 1981), xlab = "Age",
            ylab = "log10 death rate")

## -----------------------------------------------------------------------------
blocks = split_curves(x, train = 0.6, valid = 0.2)
sapply(blocks, ncol)

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
fit$d
fit$order
names(fit)
head(fit$order_selection$table)
head(fit$calibration)

## -----------------------------------------------------------------------------
fit_ac = farfima(blocks$train, blocks$valid, band = "AC", trim = 20, seed = 1, verbose = FALSE)
fit_ac

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
fit_point = farfima(blocks$train, d = "lw", trim = 20, verbose = FALSE)
fit_point

## ----error = TRUE-------------------------------------------------------------
try({
farfima(blocks$train, order = "cv", trim = 20, verbose = FALSE)
})

## -----------------------------------------------------------------------------
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)

## ----fig.height = 4-----------------------------------------------------------
plot(fit, xlab = "Age", ylab = "log10 death rate")

## ----fig.height = 4-----------------------------------------------------------
plot(fit, type = "residuals", xlab = "Age")

## -----------------------------------------------------------------------------
fc_future = predict(fit)
fc_future

## -----------------------------------------------------------------------------
fc = predict(fit, test = blocks$test)
fc

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
fc_ac5 = predict(fit, test = blocks$test, h = 5, band = "AC", verbose = FALSE)
fc_ac5$evaluation$interval

## -----------------------------------------------------------------------------
fc_point = predict(fit_point, test = cbind(blocks$valid, blocks$test), h = 5)
fc_point

## -----------------------------------------------------------------------------
fc_wrong = predict(fit_point, test = blocks$test, h = 5, verbose = FALSE)

## -----------------------------------------------------------------------------
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), ]

## ----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")

## ----fig.height = 4-----------------------------------------------------------
plot(fc, h = 1, band = "AC", xlab = "Age", ylab = "log10 death rate")

## ----fig.height = 4-----------------------------------------------------------
plot(fc, h = 3, target_test_curve = 11, band = "AC", level = 0.95, xlab = "Age", ylab = "log10 death rate")

## ----fig.height = 7-----------------------------------------------------------
plot(fc, type = "horizons", xlab = "Age", ylab = "log10 death rate")

## ----fig.height = 3.5---------------------------------------------------------
plot(fc, type = "series", at = c(0, 20, 50, 80), band = "AC", main = "Age", xlab = "Year", ylab = "log10 death rate")

## ----fig.height = 4.5---------------------------------------------------------
plot(fc, type = "recent", xlab = "Age", ylab = "log10 death rate")

## -----------------------------------------------------------------------------
head(as.data.frame(fc), 3)
df_test = as.data.frame(fc, test = TRUE)
head(df_test, 3)

## -----------------------------------------------------------------------------
ac95 = subset(df_test, band == "AC" & level == 0.95 & h == 1)
mean(ac95$observed >= ac95$lower & ac95$observed <= ac95$upper)

## -----------------------------------------------------------------------------
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)

## -----------------------------------------------------------------------------
y = frac_diff(blocks$train - rowMeans(blocks$train), d = m$d)[, -(1:20)]
dim(y)

## -----------------------------------------------------------------------------
z = frac_diff(blocks$train, d = m$d)
max(abs(frac_diff(z, d = -m$d) - blocks$train))

## -----------------------------------------------------------------------------
o_ic = select_order(y, method = "ic_full")
o_ic
select_order(y, method = "bj", seed = 1)

## -----------------------------------------------------------------------------
o_auto = select_order(y, method = "ic_auto")
select_order(y, method = "ic_auto", p = 0:3, q = 0:3)

## -----------------------------------------------------------------------------
fit_cv = farfima(blocks$train, blocks$valid, d = m$d, order = "cv", h = 5, trim = 20, verbose = FALSE)
fit_cv$order

## ----fig.width = 9, fig.height = 4--------------------------------------------
oldpar = par(mfrow = c(1, 2))
plot(o_ic)
plot(fit_cv$order_selection)
par(oldpar)

## -----------------------------------------------------------------------------
op = estimate_operators(y, order = o_ic)
op
all.equal(op[c("phi", "theta")], fit[c("phi", "theta")])

## ----fig.width = 6, fig.height = 4.5------------------------------------------
plot(op, xlab = "Age v", ylab = "Age u")

## -----------------------------------------------------------------------------
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

## -----------------------------------------------------------------------------
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)