//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  Hao Peng, John Miller
 *  @version 1.6
 *  @date    Sat Jun 13 01:27:00 EST 2017
 *  @see     LICENSE (MIT style license file).
 *
 *  @see     http://en.wikipedia.org/wiki/Autoregressive%E2%80%93moving-average_model
 *  @see     http://www.emu.edu.tr/mbalcilar/teaching2007/econ604/lecture_notes.htm
 *  @see     http://www.stat.berkeley.edu/~bartlett/courses/153-fall2010
 *  @see     http://www.princeton.edu/~apapanic/ORFE_405,_Financial_Time_Series_%28Fall_2011%29_files/slides12-13.pdf
 *  @see     www.jstatsoft.org/article/view/v027i03/v27i03.pdf
 *
 *  @title   Auto-Regressive, Moving Average (ARMA) Model
 */

package scalation.analytics
package forecaster

import scala.collection.mutable.Set
import scala.math.{log, max, min, sqrt}

import scalation.linalgebra._
import scalation.math.{double_exp, noDouble}
import scalation.minima.QuasiNewton
import scalation.plot.Plot
import scalation.random.{Normal, Uniform}
import scalation.util.banner

import Fit._
import ForecasterVec.{MAX_LAGS, testInSamp}

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMA` companion object maintains values for hyper-parameters.
 */
object ARMA
{
    /** Base hyper-parameter specification for `ARMA` as well as `AR` and `MA`
     */
    val hp = new HyperParameter
    hp += ("p", 1, 1)                                             // for the AR part
    hp += ("q", 1, 1)                                             // for the MA part
    hp += ("mle", 3, 3)                                           // optimizer to use

} // ARMA object

import ARMA._

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMA` class provides basic time series analysis capabilities for Auto-
 *  Regressive 'AR' and Moving-Average 'MA' models.  In an 'ARMA(p, q)' model,
 *  'p' and 'q' refer to the order of the Auto-Regressive and Moving-Average
 *  components of the model.  Given time series data stored in vector 'y', its
 *  next value 'y_t = y(t)' may be predicted based on prior values of 'y' and its noise:
 *  <p>
 *      y_t = δ + Σ(φ_i y_t-i) + Σ(θ_i e_t-i) + e_t
 *  <p>
 *  where 'δ' is a constant, 'φ' is the auto-regressive coefficient vector,
 *  'θ' is the moving-average coefficient vector, and 'e' is the noise vector.
 *------------------------------------------------------------------------------
 *  @param y  the original input vector (time series data)
 */
class ARMA (y: VectoD)
      extends ForecasterVec (y, -1, null)                     // call resetDF in eval_
{
    private   val DEBUG  = true                               // debug flag
    protected var differenced = false                         // flag indicating whether differencing will be applied
    protected var mu     = -0.0                               // sample mean (-0.0 means unassigned)
    protected var μ      = -0.0                               // population mean estimated using MLE
    protected var sig2   = -0.0                               // sample variance
    protected var σ2     = -0.0                               // population variance estimated using MLE

    protected var p      = 0                                  // AR order
    protected var q      = 0                                  // MA order
    protected var cap    = 0                                  // max of p and q
    protected var params = 0                                  // number of parameters estimated

    protected var zp: VectoD = null                           // vector of centered predicted/fitted values

    protected var φ:  VectoD = null                           // AR(p) coefficients
    protected var θ:  VectoD = null                           // MA(q) coefficients
    protected var δ:  Double = -0.0                           // drift/constant term

    init (y)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Initialize variables based on working time-series 'v'.
     *  Note, 'z', the centered time series, is computed later based on the MLE estimated mean μ.
     *  @param v  the working vector/time-series
     */
    def init (v: VectoD)
    {
        zp   = new VectorD (v.dim)                            // predicted/fitted values prior to uncentering
        e    = new VectorD (v.dim)                            // vector of errors/residuals
        mu   = v.mean                                         // sample mean
        sig2 = v.variance                                     // sample variance
    } // init

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the model name including its current hyper-parameters, e.g., ARMA(2, 1).
     */
    override def modelName: String = s"ARMA($p, $q)"

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector (φ concatenated with θ).
     */
    override def parameter: VectoD = φ ++ θ

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Show estimates for parameters.
     */
    def showParameterEstimates ()
    {
        println (s"differenced = $differenced")
        println (s"φ    = $φ")                                // AR parameters
        println (s"θ    = $θ")                                // MA parameters
        println (s"δ    = $δ")                                // drift
        println (s"mu   = $mu")                               // sample mean
        println (s"μ    = $μ")                                // MLE mean
        println (s"sig2 = $sig2")                             // sample variance
        println (s"σ2   = $σ2")                               // MLE variance
    } // showParameterEstimates

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Set values for the model orders 'p' and 'q'.
     *  Note, intercept/mean counts as a parameter for the time series.
     *  @param pq  the vector of model orders
     */
    def setPQ (pq: VectoI): ARMA =
    {
        val n = pq.dim
        if (n > 0) p = pq(0)
        φ = new VectorD (p)                                   // AR coefficients
        if (n > 1) q = pq(1)
        θ = new VectorD (q)                                   // MA coefficients
        cap    = max (p, q)                                   // greatest lag
        params = p + q + (if (differenced) 0 else 1)          // number of parameters
        this
    } // setPQ

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train/fit an `ARMA` model to the times series data.  Must call 'setPQ' first.
     *  It uses BFGS, a Quasi-Newton optimizer, to minimize the negative log-likelihood.
     */
    def train (): ARMA =
    {
        super.train (null, y)
        val solver  = new QuasiNewton (nll)                   // nonlinear optimizer
        val b       = new VectorD (params + 1)                // parameter values incl. σ2

        if (! differenced) b(b.size-2) = mu
        b(b.size-1) = sqrt (sig2)                             // sample standard deviation, initial est. for σ parameter
        solver.solve (b)                                      // find b that minimizes negative log-likelihood

        δ = μ * (1 - φ.sum)                                   // update drift value
//      δ = stats.mu * (1 - φ.sum)

        if (DEBUG) showParameterEstimates ()
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** The negative log-likelihood function to be minimized.
     *  @see math.unice.fr/~frapetti/CorsoP/Chapitre_4_IMEA_1.pdf, page 36
     *  @see spia.uga.edu/faculty_pages/monogan/teaching/ts/Barima.pdf
     *  @see stats.stackexchange.com/questions/77663/arima-estimation-by-hand
     *  @param b  the input parameter vector
     */
    def nll (b: VectoD): Double =
    {
        if (b.size != params + 1) flaw ("nll", "input parameter vector size incorrect")
        for (i <- 0 until p)   φ(i)   = b(i)
        for (i <- p until p+q) θ(i-p) = b(i)
        if (! differenced) μ = b(b.size-2)
        σ2 = b.last~^2

        updateFittedValues ()
    } // nll

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Update the vector of fitted values 'zp', the vector of errors 'e', and
     *  return the negative log-likelihood '-ll'.
     *  @see `Fit` for definition of 'll'.
     */
    def updateFittedValues (): Double =
    {
        if (! differenced) for (i <- z.range) z(i) = y(i) - μ    // for undifferenced time series, center using est. μ

        zp(0) = z(0)                                             // no past values or errors => copy actual
        for (t <- 1 until zp.dim) {
            e(t-1)  = z(t-1) - zp(t-1)                           // error in previous forecast
            var sum = 0.0
            for (j <- φ.range if t-j > 0) sum += φ(j) * z(t-1-j)
            for (j <- θ.range if t-j > 0) sum += θ(j) * e(t-1-j)
            zp(t) = sum
        } // for

        -ll (e.normSq / m, σ2, m)                                // return negative log likelihood
    } // updateFittedValues

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the vector of predicted/fitted values on the training/full dataset.
     *  Based on 'zp' calculated in the 'updateFittedValues' method.
     */
    override def predictAll (): VectoD = zp + μ

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the vector of zero-centered predicted/fitted values on the training/full
     *  dataset.  Based on 'zp' calculated in the 'updateFittedValues' method.
     */
    override def predictAllz (): VectoD = predictAll () - μ

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Produce h-steps-ahead forecast for ARMA models.
     *  @see ams.sunysb.edu/~zhu/ams586/Forecasting.pdf
     *  @param t  the time point from which to make forecasts (in the original scale)
     *  @param h  the number of steps to forecast, must be at least one.
     */
    override def forecast (t: Int = y.dim, h: Int = 1): VectoD =
    {   
        if (t > y.dim) flaw ("forecast", s"t ($t) cannot be greater than y.dim (${y.dim})") 
        
        val zf = new VectorD (cap + h)                           // forecasted centered values
        val e_ = new VectorD (cap + h)                           // available observed errors
        
        for (i <- 0 until cap if t-cap+i >= 0) {                 // seed with first cap = max(p, q) values
            zf(i) = z(t-cap+i)                                   // copy first cap values
            e_(i) = e(t-cap+i)                                   // unveil first cap errors (observed in training)
        } // for

        for (i <- cap until zf.dim) {                            // start at t = cap (enough for first value to forecast)
            var sum = 0.0
            for (j <- φ.range) sum += φ(j) * zf(i-1-j)
            for (j <- θ.range) sum += θ(j) * e_(i-1-j)
            zf(i) = sum
        } // for
        val f = zf.slice (cap)                                   // dump first cap forecasts
        f + μ                                                    // return the vector of forecasts
    } // forecast

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Forecast values for all time points using 1 through 'h'-steps ahead forecasts.
     *  The 'h'-th row of matrix is the horizon 'h' forecast (where 'h = 0' is actual data).
     *  @param h  the forecasting horizon, number of steps ahead to produce forecasts
     */
    override def forecastAll (h: Int = 1): MatriD =
    {   
        val yf  = new MatrixD (y.dim, h+1)                       // forecasts for all horizons h & time points t
        yf.setCol(0, y)                                          // first row is actual values
        val cut = cap                                            // cut over from actual to forecasted values
        
        for (t <- y.range) {
            if (t < cut) {
                for (k <- 1 to h) yf(t, k) = y(t)                // copy first cut observed values from y
            } else {
                val ft = forecast (t, h)                         // forecasts at time point t, horizons 1 to h
                for (k <- 1 to h if t+k-1 < y.dim) {             
                    yf(t+k-1, k) = ft(k-1)                       // place forecasts diagonally
                } // for
            } // if
        } // for
        
        // fill in blank values in first few rows where no forecasts can be produced by copying values from previous columns
        for (k <- 2 to h; t <- cut until cut+k-1) yf(t, k) = yf(t, k-1)   // copy forecasted values
        
        yf                                                       // return matrix of forecasted values
    } // forecastAll

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Forecast values for all 'm' time points and all horizons (1 through 'h'-steps ahead).
     *  Record these in the 'yf' matrix, where
     *  <p>
     *      yf(t, k) = k-steps ahead forecast for y_t
     *  <p>
     *  Note, 'yf.col(0)' is set to 'y' (the actual time-series values).
     *  Do not forecast errors, rather use observed errors from training and make sure not
     *  to use errors that would correspond to knowing future errors (all future errors should
     *  be assumed to be 0).
     *  @see https://otexts.com/fpp3/arima-forecasting.html, section 9.8
     *  @param h  the maximum forecasting horizon, number of steps ahead to produce forecasts
     */
    def forecastAll2 (h: Int): MatriD =
    {
        val zf = new MatrixD (z.dim, h+1)                        // forecast matrix: rows - time, cols - horizon
        zf.setCol (0, z)                                         // first column is actual values, horizon 0

        for (k <- 1 to h) {                                      // loop through k-steps ahead forecasts
            val e_ = new VectorD (z.dim)                         // redetermine errors from a clean slate

            for (t <- 0 until cap) {                             // seed the first cap = max(p, q) values
                zf(t, k) = z(t)                                  // copy first cap actual values
                e_(t)    = e(t)                                  // copy first cap errors (observed in training)
            } // for

            for (t <- cap until z.dim) {                         // forecast from cap to the end of time-series
                if (t-k >= 0) e_(t-k) = e(t-k)                   // unveil previous error at time t-k
                var sum = 0.0
                for (j <- φ.range if t-j > 0) sum += φ(j) * zf(t-1-j, max (0, k-1-j))
                for (j <- θ.range if t-j > 0) sum += θ(j) * e_(t-1-j)
                zf(t, k) = sum                                   // centered forecast for time t
            } // for
        } // for
        zf + μ                                                   // return uncentered forecasts
    } // forecastAll2

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Obtain residuals/errors in the original scale.
     */
    def residuals: VectoD = e

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Compute the error (difference between actual and predicted) and useful
     *  diagnostics for the dataset.
     *  @param e    vector of residuals
     *  @param y    vector of observed values
     *  @param yp   vector of predicted values
     */
    def eval_ (e: VectoD = residuals, y: VectoD = y, yp: VectoD = predictAll ()): ARMA =
    {
        resetDF (params, y.dim - params)
        diagnose (e, y, yp)
        this
    } // eval_

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Perform forward selection to find the most predictive variables to add the
     *  existing model, returning the variables to add and the new model.
     *  Note, all lags up and including 'p' define the model.
     *  FIX - select subsets of the lags, e.g., Set (1, 2, 5)
     *  @see `Fit` for index of QoF measures.
     *  @param cols     the lags/columns currently included in the existing model (ignored)
     *  @param index_q  index of Quality of Fit (QoF) to use for comparing quality
     */
    def forwardSel (cols: Set [Int] = null, index_q: Int = index_rSq): (Int, ARMA) =
    {
        val rSq  = new MatrixD (MAX_LAGS, 3)                  // R^2, R^2 Bar, R^2 cv
        var p_mx   = -1                                       // best column, so far
        var mod_mx = null.asInstanceOf [ARMA]                 // best model, so far
        var fit_mx = noDouble                                 // best fit, so far

        for (p <- 1 to MAX_LAGS) {
            setPQ (VectorI (p, p)).train ().eval_ ()          // set p and q, retrain model, evaluate QoF
            val fit_p = fit(index_q)                          // new fit
            if (fit_p > fit_mx) { p_mx = p; mod_mx = this; fit_mx = fit_p }
            rSq(p-1)  = VectorD (fit_p, -1.0, -1.0)           // use new model, mod_p, no cross
        } // for
        if (p_mx == -1) {
            flaw ("forwardSel", "could not add additional lags")
        } // if
        println (s"forwardSel: add $p_mx giving qof = $fit_mx")
        println (s"rSq = $rSq")
        new Plot (null, rSq.col(0), null, s"forwardSel $modelName", lines = true)
        (p_mx, mod_mx)
    } // forwardSel

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Perform forward selection to find the most predictive variables to add the
     *  existing model, returning the variables to add and the new model.
     *  Note, all lags up and including 'p' define the model.
     *  FIX - select subsets of the lags, e.g., Set (1, 2, 5)
     *  @see `Fit` for index of QoF measures.
     *  @param step     the amount to increment p and q for each iteration
     *  @param cols     the lags/columns currently included in the existing model (ignored)
     *  @param index_q  index of Quality of Fit (QoF) to use for comparing quality
     */
    def forwardSel2 (step: VectoI = VectorI (1, 1), cols: Set [Int] = null,
                     index_q: Int = index_rSq): (Int, ARMA) =
    {
        val rSq    = new MatrixD (MAX_LAGS, 3)                // R^2, R^2 Bar, R^2 cv
        var p_mx   = -1                                       // best column, so far
        var mod_mx = null.asInstanceOf [ARMA]                 // best model, so far
        var fit_mx = noDouble                                 // best fit, so far

        for (p <- 1 to MAX_LAGS) {
            val pq = step * p
            setPQ (pq).train ().eval_ ()                      // set p and q, retrain model, evaluate QoF
            val fit_p = fit(index_q)                          // new fit
            if (fit_p > fit_mx) { p_mx = p; mod_mx = this; fit_mx = fit_p }
            rSq(p-1)  = VectorD (fit_p, -1.0, -1.0)           // use new model, mod_p, no cross
        } // for
        if (p_mx == -1) {
            flaw ("forwardSel2", "could not add additional lags")
        } // if
        println (s"forwardSel2: add $p_mx giving qof = $fit_mx")
        println (s"for step = $step, rSq = $rSq")
        new Plot (null, rSq.col(0), null, s"forwardSel2 $modelName, step = $step", lines = true)
        (p_mx, mod_mx)
    } // forwardSel2

} // ARMA class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMATest` object is used to test the `ARMA` class.
 *  > runMain scalation.analytics.forecaster.ARMATest
 */
object ARMATest extends App
{
    println ("ARMA!")
    val m = 100
    val noise = Uniform (-5, 5)
    val y = VectorD (for (i <- 0 until m) yield i + noise.gen)

    val arma = new ARMA (y)                          // time series data: y vs. t

    // Build ARMA(p, q) models

    arma.setPQ (VectorI (1, 0)).train ().eval_ ()
    testInSamp (y, 3, arma, 1)

    arma.setPQ (VectorI (2, 0)).train ().eval_ ()
    testInSamp (y, 3, arma, 1)

    arma.setPQ (VectorI (0, 3)).train ().eval_ ()
    testInSamp (y, 3, arma, 1)

    println (arma.forecastAll (3))

    arma.plotFunc2 (arma.acF, "ACF")                 // must call super.train, so that ACF is actually computed - FIX
    arma.plotFunc2 (arma.pacF, "PACF")               // must call super.train, so that PACF is actually computed

} // ARMATest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMATest2` object is used to test the `ARMA` class.
 *  > runMain scalation.analytics.forecaster.ARMATest2
 */
object ARMATest2 extends App
{
    val m = 20
    val noise = Normal (0, 2)
    val y = VectorD (for (i <- 0 until m) yield i + noise.gen)

    println (s"y = $y")

    val p = 1
    val q = 1
    val steps = 2                                    // number of steps for the forecasts

    val arma = new ARMA (y)                          // time series data: y vs. t

    // Build ARMA(p, q) models

    arma.setPQ (VectorI (p, 0)).train ().eval_ ()
    testInSamp (y, p, arma, 1)
    var yf = arma.forecast (t = m, h = steps)
    println (s"$steps-step ahead forecasts using ARMA($p, 0) model = $yf")

    arma.setPQ (VectorI (0, q)).train ().eval_ ()
    testInSamp (y, q, arma, 1)
    yf = arma.forecast (t = m, h = steps)
    println (s"$steps-step ahead forecasts using ARMA(0, $q) model = $yf")

    arma.setPQ (VectorI (p, q)).train ().eval_ ()
    testInSamp (y, p+q, arma, 1)
    yf = arma.forecast (t = m, h = steps)
    println (s"$steps-step ahead forecasts using ARMA($p, $q) model = $yf")

} // ARMATest2 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMATest3` object is used to test the `ARMA` class on real data:
 *  Forecasting lake levels.
 *  @see cran.r-project.org/web/packages/fpp/fpp.pdf
 *  > runMain scalation.analytics.forecaster.ARMATest3
 */
object ARMATest3 extends App
{
    import ForecasterVec.y

    for (h <- 1 to 2) {                                        // forecasting horizon
        for (p <- 0 to 4) {                                    // auto-regressive hyper-parameter settings
            for (q <- 0 to 4 if p+q > 0) {                     // moving-average hyper-parameter settings
                banner (s"Test3: ARMA ($p, $q) with h = $h")
                val mod = new ARMA (y)                         // create an ARMA model
                mod.setPQ (VectorI (p, q)).train ().eval_ ()   // set p and q, train and evaluate the ARMA model
                println (mod.report)
                val yp  = mod.predictAll ()                    // predicted/fitted values
                val yfa = mod.forecastAll (h)
                val yfb = mod.forecastAll2 (h)
                assert (yfa == yfb)
                val yf  = yfa.col(h)                           // forecasted values - h steps ahead
                val yf2 = yfb.col(h)                           // forecasted 2 values - h steps ahead

                new Plot (null, y, yp, s"Plot of y & yp, predicted ARMA ($p, $q) vs. t", true)
                new Plot (null, y, yf, s"Plot of y & yf, forecasted (h = $h) ARMA ($p, $q) vs. t", true)
                new Plot (null, y, yf2, s"Plot of y & yf2, forecasted2 (h = $h) ARMA ($p, $q) vs. t", true)

                if (h == 1) ForecasterVec.differ (yp, yf, allow = true)
                ForecasterVec.differ (yf, yf2, allow = true)
                val skip = max (p, q)                          // skip the cap start-up

                banner (s"ARMATest3: QoF (@h = $h) for yf = mod.forecastAll")
                println (s"rSq (yf)   for h = $h is ${rSqF (y, yf)}")
                println (s"rSq (yf)   for h = $h is ${rSqF (y, yf, max (p, q))}, skip = $skip")
                println (s"sMAPE (yf) for h = $h is ${smapeF (y, yf)}")
                println (s"sMAPE (yf) for h = $h is ${smapeF (y, yf, max (p, q))}, skip = $skip")

                banner (s"ARMATest3: QoF (@h = $h) for yf2 = mod.forecastAll2")
                println (s"rSq (yf2)   for h = $h is ${rSqF (y, yf2)}")
                println (s"rSq (yf2)   for h = $h is ${rSqF (y, yf2, max (p, q))}, skip = $skip")
                println (s"sMAPE (yf2) for h = $h is ${smapeF (y, yf2)}")
                println (s"sMAPE (yf2) for h = $h is ${smapeF (y, yf2, max (p, q))}, skip = $skip")
            } // for
        } // for
    } // for

} // ARMATest3 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMATest4` object is used to test the `ARMA` class.
 * > runMain scalation.analytics.forecaster.ARMATest4
 */
object ARMATest4 extends App
{
    val path = BASE_DIR + "travelTime.csv"
    val data = MatrixD (path)

    val y = data.col(1)

    val p = 1
    val q = 1
    val steps = 1                                  // number of steps for the forecasts

    val arma = new ARMA (y)                        // time series data: y vs. t

    println (s"y = $y")

    // Build ARMA(p, q) models

    arma.setPQ (VectorI (p, 0)).train ().eval_ ()
    testInSamp (y, p, arma, 1)
    var yf = arma.forecast (t = y.dim, h = steps)
    println (s"$steps-step ahead forecasts using ARMA($p, 0) model = $yf")

    arma.setPQ (VectorI (0, q)).train ().eval_ ()
    testInSamp (y, q, arma, 1)
    yf = arma.forecast (t = y.dim, h = steps)
    println (s"$steps-step ahead forecasts using ARMA(0, $q) model = $yf")

    arma.setPQ (VectorI (p, q)).train ().eval_ ()
    testInSamp (y, p+q, arma, 1)
    yf = arma.forecast (t = y.dim, h = steps)
    println (s"$steps-step ahead forecasts using ARMA($p, $q) model = $yf")

} // ARMATest4 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMATest5` object is used to test the `ARMA` class.
 *  > runMain scalation.analytics.forecaster.ARMATest5
 */
object ARMATest5 extends App
{
    val m     = 50
    val steps = 2                                                  // number of steps for the forecasts
    val sig2  = 10000.0
    val noise = Normal (0.0, sig2)
    val y = VectorD (for (i <- 0 until m) yield 40 * (i-1) - (i-2) * (i-2) + noise.gen)
    
    val arma = new ARMA (y)                                        // time series model ARMA
    
    banner ("Build ARMA(1, 0) model")                           
    arma.setPQ (VectorI (1, 0)).train ().eval_ ()                  // train for ARMA(1, 0) model
    testInSamp (y, 2, arma, 1)
    
    for (p <- 1 to 3) { 
        banner (s"Build ARMA($p, $p) model")
        arma.setPQ (VectorI (p, p)).train ().eval_ ()              // retrain for ARMA(p, p) model
        testInSamp (y, p+p, arma, 1)
        
        banner ("Make Forecasts")
        val yf = arma.forecast (steps)
        println (s"$steps-step ahead forecasts using ARMA($p, $p) model = $yf")
    } // for

} // ARMATest5 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARMATest6` object is used to test the `ARMA` class on real data:
 *  Forecasting lake levels.  Select the best number of lags.
 *  @see cran.r-project.org/web/packages/fpp/fpp.pdf
 *  > runMain scalation.analytics.forecaster.ARMATest6
 */
object ARMATest6 extends App
{
    import ForecasterVec.y

    val arma = new ARMA (y)                                        // create model for time series data
    arma.setPQ (VectorI (1, 1)).train ().eval_ ()                  // train the model and evaluate
    val res = arma.forwardSel ()
    println (s"forwardSel: $res")

    for ((sp, sq) <- Array ((1, 0), (2, 1), (1, 1), (1, 2), (0, 1))) {
        val res = arma.forwardSel2 (VectorI (sp, sq))
        println (s"forwardSel2: $res")
    } // for

} // ARMATest6 object