\documentclass{article}
\usepackage[margin=1in]{geometry}
\usepackage{hyperref}
\usepackage{Sweave}
%\VignetteIndexEntry{Introduction to Local Latent Flow Forecasting}
%\VignetteEngine{utils::Sweave}
%\VignetteEncoding{UTF-8}
\title{Forecasting with flowz}
\author{Giancarlo Vercellino}
\date{}
\begin{document}
\SweaveOpts{concordance=TRUE}
\maketitle
<<setup, echo=FALSE, results=hide>>=
flowz_old_options <- options(width = 72)
library(flowz)
@
\section{Inputs and a complete first forecast}
The package models a joint distribution of future levels of a regularly
observed univariate series. Supply a finite numeric vector, ordered from
oldest to newest observation. A univariate \texttt{ts} is accepted, but its
calendar attributes are not carried into the output. One forecast horizon
means one observation interval. Resolve missing data before fitting; dropping
rows can change the meaning of that interval.

The following deterministic series needs no external data. These small
computational settings make the vignette fast to rebuild; they are not
accuracy recommendations for a production application.
<<data>>=
t <- seq_len(260)
x <- cumsum(sin(t / 4) + cos(t / 11))
@
<<fit>>=
fit <- flowz(
  x, horizon = 3, transform = "none",
  baseline = "ewma", baseline_control = list(), autotune = TRUE,
  drift_lambda = 0.90, vol_lambda = 0.94,
  min_train = NULL, innovation_lags = 8,
  n_hidden = 10, lambda = 1e-3, flow_aug = 1,
  n_steps = 3, n_sims = 30,
  cal_fraction = 0.20, test_fraction = 0.20,
  cal_sims = 30, test_sims = 30,
  spectral_modes = NULL, spectral_spar = 0.65,
  spectral_shrink_to_one = 0.10,
  spectral_tau_bounds = c(0.15, 2.50),
  distribution_bw_factor = 0.50,
  seed = 123, n_origins = NULL,
  parallel = FALSE, n_cores = NULL,
  plot = FALSE, verbose = FALSE
)
fit$forecast$calibrated$summary
dim(fit$forecast$calibrated$samples)
@
Each row of the sample matrix is one simulated future path; columns represent
horizons. The primary model is \texttt{calibrated}, after spectral calibration.
The \texttt{uncalibrated} model contains paths before calibration, and
\texttt{baseline} contains the selected baseline forecast. Summary quantiles describe the
empirical simulated paths. They can differ slightly from quantiles returned
by the smooth marginal distribution functions.

\section{A consistent output structure}
The fitted object has five visible sections: \texttt{forecast},
\texttt{evaluation}, \texttt{calibration}, \texttt{model}, and
\texttt{metadata}. The forecast contains three model components:
\texttt{calibrated}, \texttt{baseline}, and \texttt{uncalibrated}.
Every model component uses the same structure:
\begin{description}
\item[summary] One data-frame row per horizon; means, medians,
standard deviations, and quantiles are on the original data scale.
\item[samples] A path matrix whose columns are named \texttt{h1},
\texttt{h2}, and so on.
\item[pred\_fun] Marginal functions grouped as
\texttt{dfun}, \texttt{pfun}, \texttt{qfun}, and \texttt{rfun} horizon lists.
\end{description}
<<output-structure>>=
fc <- fit$forecast
names(fc)
fc$calibrated$pred_fun$qfun$h3(c(.05, .5, .95))
fc$baseline$pred_fun$qfun$h3(c(.05, .5, .95))
fc$uncalibrated$pred_fun$qfun$h3(c(.05, .5, .95))
fit$evaluation$metadata$scale
@
The \texttt{calibrated} group is the final calibrated forecast.
The \texttt{baseline} group is the selected EWMA, ARIMA, ARFIMA, or structural
model, and \texttt{uncalibrated} contains Flow marginals before spectral
correction. Each model owns all four
functions: \texttt{dfun}, \texttt{pfun}, \texttt{qfun}, and \texttt{rfun}.
Evaluation belongs to the fitted model, not to an updated forecast whose
future observations are unknown. Use \texttt{fit\$evaluation} for scores,
\texttt{fit\$calibration} for spectral estimates, and
\texttt{fit\$metadata\$diagnostics} for fitting details. There are no
duplicated forecast aliases: select a model first, then its summary, paths,
or marginal functions. See \texttt{?flowz\_forecast} for the exact contract.

\section{Choosing a baseline engine}
The baseline supplies the causal mean trajectory and horizon-specific
uncertainty used to define Flow residuals. Four engines are available:
\begin{description}
\item[ewma] Local exponentially weighted drift and volatility; the default.
\item[arima] Nonseasonal ARIMA. Autotuning selects the order on the initial
training window.
\item[arfima] Fractionally integrated ARMA for persistent, long-memory series.
\item[structural] A structural local-level or local-trend state-space model.
\end{description}
The selected specification is refitted using only observations available at
each rolling origin. This prevents future observations from entering baseline
estimation. EWMA states remain conditioning features for the Flow layer across
all four choices, so only the baseline trajectory and scale change.

The following calls illustrate the interface; they are not run during the
vignette build because rolling ARFIMA estimation is more costly.
<<baseline-options, eval=FALSE>>=
fit_arima <- flowz(x, baseline = "arima", autotune = TRUE)
fit_arfima <- flowz(x, baseline = "arfima", autotune = TRUE)
fit_structural <- flowz(x, baseline = "structural",
  baseline_control = list(type = "trend"), autotune = FALSE)
fit_arima$metadata$baseline_label
@
For a fixed ARIMA specification, use for example
\texttt{baseline\_control = list(order = c(1,1,1))} with
\texttt{autotune = FALSE}. ARFIMA accepts \texttt{drange} and
\texttt{estim}; the structural engine accepts \texttt{type = "level"} or
\texttt{"trend"}. An isolated fitting failure is recorded in
\texttt{fit\$metadata\$diagnostics\$baseline} and falls back to Drift/EWMA.

\section{Method and interpretation}
Write $y_t=x_t$ for an untransformed fit and $y_t=\log(x_t)$ for a log fit.
For the default EWMA engine, let $\Delta_t=y_t-y_{t-1}$. It tracks local
drift and innovation variance using
\[
 \mu_t=\lambda_\mu\mu_{t-1}+(1-\lambda_\mu)\Delta_t,\qquad
 e_t=\Delta_t-\mu_{t-1},
\]
\[
 \sigma_t^2=\lambda_\sigma\sigma_{t-1}^2+
                  (1-\lambda_\sigma)e_t^2.
\]
At origin $i$, its horizon-$h$ baseline is $y_i+h\mu_i$, with scale
$\sigma_i\sqrt{h}$. The other engines replace this trajectory and scale with
their causal forecasts. Horizon corrections estimated from training origins
standardize the multi-step level errors for every engine.

The model learns the joint error vector directly. Conditioning features
include lagged standardized innovations, drift, volatility, and momentum.
Random Gaussian starting points and training error vectors are interpolated
linearly at random times. An extreme learning machine (ELM), with random
tanh hidden features and ridge-fitted output weights, approximates the
conditional velocity field. Fourth-order Runge--Kutta integration generates
future residual paths, which are added to the baseline.

Calibration uses a discrete cosine transform (DCT) of path deviations around
the Flow sample mean. Calibration-origin errors and simulated variation
determine spectral dispersion multipliers. Smoothing, shrinkage towards one,
and bounds regularize those multipliers. Unselected higher-frequency modes
are left unchanged before the final bounds are applied.
Calibration preserves the sample mean on the fitting scale; nonlinear
back-transformation need not preserve the mean on the price scale.

There is a single fitted velocity field. The simulated uncertainty represents
the learned residual distribution, without a bootstrap ensemble for parameter
uncertainty. Coverage is empirical, not an exact finite-sample guarantee.

\section{Evaluation and its limits}
Eligible historical origins begin at \texttt{min\_train} and end at
\texttt{length(x) - horizon}. Automatic origin selection uses a heuristic
with a cap of 500; \texttt{n\_origins = Inf} requests all eligible origins.
Selected origins are split chronologically into training, calibration, and
test blocks. Calibration and test blocks each contain at least ten origins;
at least twenty must remain for training. A small manually specified origin
count can therefore fail even when the input series is long.

The current implementation does not enforce an embargo between forecast
windows. Consecutive origins, including origins on opposite sides of a split,
may have overlapping future observations. Chronological ordering alone does
not make these scores independent holdout estimates. Inspect the spacing and
reserve a separate future period when assessing generalization.
<<evaluation>>=
o <- fit$metadata$diagnostics$origin_selection
c(available = o$available_origins,
  selected = o$selected_origins,
  median_spacing = o$median_spacing)
sc <- fit$calibration
c(training = sc$train_origins,
  calibration = sc$calibration_origins,
  test = sc$test_origins)
fit$evaluation$metrics$aggregate
@
RMSE and MAE measure errors of point forecasts. CRPS assesses marginal
distribution accuracy. Smaller values are preferable for all three.
Empirical 90 percent coverage measures the fraction of observations inside
the 5th--95th percentile interval; coverage alone does not measure how useful
or narrow that interval is. Compare Flow against the selected baseline across horizons,
not only by an aggregate score.
<<horizon-scores>>=
head(fit$evaluation$metrics$per_horizon)
@
Changing settings repeatedly based on these test scores makes this block a
model-selection set. Use a further untouched future period for final
assessment. Serial dependence and overlapping horizons also affect the
interpretation of pooled calibration diagnostics.

\section{Marginal probabilities and joint events}
The named horizon lists make marginal queries explicit:
<<marginals>>=
fun <- fit$forecast$calibrated$pred_fun
fun$qfun$h3(c(0.05, 0.50, 0.95))
threshold <- fun$qfun$h3(0.5)
fun$pfun$h3(threshold)
fun$dfun$h3(threshold)
set.seed(42)
head(fun$rfun$h3(10))
@
For horizon $h$, Flow uses a Gaussian-kernel mixture centered on simulated
values. The bandwidth multiplier controls smoothing of these marginal
functions, without altering stored paths or spectral calibration. The
\texttt{dfun}, \texttt{pfun}, and \texttt{qfun} functions accept the usual
\texttt{log}, \texttt{lower.tail}, and \texttt{log.p} arguments where applicable.

Use complete sample rows for joint events. Independent calls to
\texttt{rfun\$h1} and \texttt{rfun\$h2} do not retain cross-horizon dependence.
For example, estimate the probability that all three future levels exceed
the last observed level:
<<joint-event>>=
paths <- fit$forecast$calibrated$samples
mean(rowSums(paths > tail(x, 1)) == ncol(paths))
@
With only 30 simulated paths this is a coarse Monte Carlo estimate.
Increase the simulation budget for decisions that depend on tail probabilities.

\section{Updating the forecast}
Pass the full updated history. Prediction recomputes the local state and
context, keeping the fitted horizon, velocity field, and calibration fixed.
It does not retrain the model.
<<prediction>>=
t_new <- seq_len(265)
x_new <- cumsum(sin(t_new / 4) + cos(t_new / 11))
pr <- predict(fit, x_new, n_sims = 50, seed = 42)
pr$forecast$calibrated$summary
pr$forecast$calibrated$pred_fun$qfun$h3(c(0.05, 0.95))
@
The generic \texttt{predict()} method defaults to 1000 paths. Specify
\texttt{n\_sims} explicitly for comparable runs.
Refit when the learned residual dynamics or calibration may have changed.

\section{Positive prices and log transformation}
Use log transformation for strictly positive data when additive errors on the
log scale are appropriate. This produces positive forecast paths and marginal
quantiles under ordinary floating-point ranges.
<<prices>>=
prices <- 100 * exp(x / 20)
price_fit <- flowz(prices, horizon = 3, transform = "log",
  n_hidden = 10, flow_aug = 1, n_steps = 3,
  n_sims = 30, cal_sims = 30, test_sims = 30, parallel = FALSE)
price_fit$forecast$calibrated$summary[, c("horizon", "mean", "q05", "q95")]
price_fit$forecast$calibrated$pred_fun$qfun$h3(c(0.05, 0.50, 0.95))
stopifnot(all(price_fit$forecast$calibrated$samples > 0))
@
Forecast summary means and standard deviations are calculated from the
back-transformed samples; exponentiating a log mean does not give an
arithmetic price mean. The baseline trajectory is exponentiated directly.
Kernel bandwidths, fitted parameters, and evaluation diagnostics remain on
the log scale.
In particular, the test RMSE panel reports log-scale errors even though the
current forecast panel displays price levels.

\section{Reading the four-panel plot}
<<diagnostics, fig=TRUE, width=8, height=6, echo=TRUE>>=
plot(price_fit)
@
Calibrated Flow is thick solid black; the selected baseline is dashed gray.
The upper-left panel shows the observed history, current median forecast,
baseline trajectory, and empirical 90 percent band.
The other panels compare test RMSE, empirical coverage against the dotted
90 percent reference, and PIT quantiles against the dotted uniform reference.
The PIT panel pools horizons and test origins and may hide differences by
horizon. A line near the diagonal is useful evidence, not proof of calibration.

Set \texttt{plot = TRUE} in \texttt{flowz()} to draw this display immediately.
\texttt{plot(fit)} always uses the original fitted forecast and diagnostics.
It does not display the result of a later \texttt{predict()} call. Currently
prediction-time plotting is available for untransformed fits; for a log fit,
draw updated predictions from the returned original-scale summary or samples.

\section{Choosing computational settings}
For a reproducible first run, keep the smoothing parameters fixed and use
automatic origin selection. Increase one computational control at a time,
holding the seed and evaluation design constant. Inspect both per-horizon
scores and forecast bands after each change.
For the EWMA engine, the two persistence parameters are fixed inputs, with
defaults 0.90 and 0.94; larger values adapt more slowly. For ARIMA, ARFIMA,
and structural engines, \texttt{autotune = TRUE} selects the specification
on the initial training window. It does not tune the Flow hyperparameters.

Increase \texttt{n\_hidden} to expand capacity, \texttt{flow\_aug} to supply
more random training interpolations, and \texttt{n\_steps} to improve integration
resolution. These increase computation and do not guarantee better forecasts.
The ridge penalty regularizes the fitted output weights.
Increase \texttt{cal\_sims}, \texttt{test\_sims}, and \texttt{n\_sims} to reduce
simulation noise in their respective stages.
More spectral shrinkage gives a weaker dispersion correction; narrow bounds
limit extreme corrections. Validate these choices using an appropriate
external chronological design.

Parallel calibration and testing are enabled by default. Automatic worker
selection uses at most two cores; an explicit \texttt{n\_cores} can request
more. Windows uses PSOCK workers and other platforms use forked processes.
These vignette examples use serial execution to keep build overhead small.
For small models, process startup and communication can cost more than they
save. Use \texttt{fit\$metadata\$diagnostics\$timing} to inspect elapsed
stage times.

A fixed seed supports repeatable runs for the same data, settings, and software
environment. Fitting and prediction advance R's random-number generator.
When prediction uses \texttt{seed = NULL}, it continues from the current
generator state; provide an explicit seed to compare repeated forecasts.
Store the complete object using \texttt{saveRDS}; \texttt{readRDS} restores
the fitted object and its internal prediction state. Consult \texttt{?flowz} and
\texttt{?predict.flowz} for the full return structure and method signatures.
\section{Common questions and troubleshooting}
\begin{description}
\item[Too few origins.] The split needs at least twenty training origins
and at least ten in each later block. Start with \texttt{n\_origins = NULL}.
Longer horizons may require a longer input history.
\item[Different quantiles.] Summary quantiles describe stored paths.
The marginal \texttt{qfun} functions use a smooth kernel mixture. Small
differences are expected, especially with a small simulation budget.
\item[Changing prediction results.] Specify the same prediction seed and
simulation count for a direct comparison.
\item[Unexpected score units.] With a log fit, evaluation scores remain on
the log scale. Evaluate original-scale predictions separately when the
decision requires price-unit errors.
\item[Parallel overhead.] Serial execution can be faster for small models.
Compare stage timings with the same data and settings before increasing
the worker count.
\end{description}
<<cleanup, echo=FALSE, results=hide>>=
options(flowz_old_options)
@
\end{document}