//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  Hao Peng, John Miller, Michael Cotterell
 *  @version 1.5
 *  @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
 */

package scalation.analytics.forecaster

import scala.math.{max, min}

import scalation.analytics.{BASE_DIR, Regression}
import scalation.linalgebra.{MatriD, MatrixD, VectoD, VectorD}
import scalation.plot.Plot
import scalation.random.{Normal, Random}
import scalation.stat.vectorD2StatVector

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `AR` class provides basic time series analysis capabilities for Auto-
 *  Regressive 'AR'.  In an  'AR(p)' model, 'p' refers to the order of the
 *  Auto-Regressive components of the model.  `AR` models are often used for forecasting.
 *  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 = c + Σ(φ_i y_t-i) + Σ(θ_i e_t-i) + e_t
 *  <p>
 *  where 'c' is a constant, 'φ' is the autoregressive coefficient vector,
 *  and 'e' is the noise vector.
 *------------------------------------------------------------------------------
 *  @param t  the time vector
 *  @param y  the input vector (time series data)
 */
class AR (t: VectoD, y: VectoD)
      extends ForecasterVec (t, y)
{
    private val DEBUG = false                            // debug flag
    private val x     = (y - mu).asInstanceOf [VectorD]  // work with mean zero time series
    private val sig2  = x.variance                       // the sample variance
    private val ac    = new VectorD (ml+1)               // auto-covariance
    for (t <- ac.range) ac(t) = x acov t
            val acf   = ac / sig2                        // Auto-Correlation Function (ACF)
            var pacf: VectoD = null                      // Partial Auto-Correlation Function (PACF)

    private var p = 0                                    // order of AR component
    private var φ: VectoD = null                         // AR(p) coefficients

    if (DEBUG) {
        println ("n    = " + n)                          // number of data points
        println ("x    = " + x)                          // zero-mean time series
        println ("mu   = " + mu)                         // mean
        println ("sig2 = " + sig2)                       // variance
        println ("ac   = " + ac)                         // auto-covariance
        println ("acf  = " + acf)                        // ACF
    } // for

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Set value for 'p'.
     *  @param p_  the order of the AR part of the model
     */
    def setPQ (p_ : Int) { p = p_ }

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train/fit an `ARMA` model to times the series data.  Must call setPQ first.
     */
    def train (): AR = train (p)

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Estimate the coefficient vector 'φ' for a 'p'th order Auto-Regressive 'AR(p)'
     *  model.
     *  <p>
     *      x_t = φ_0 * x_t-1 + ... + φ_p-1 * x_t-p + e_t
     *  <p>
     *  Uses the Durbin-Levinson Algorithm to determine the coefficients.
     *  The 'φ' vector is 'p'th row of 'psi' matrix (ignoring the first (0th) column).
     *  @param p_  the order of the AR model
     */
    def train (p_ : Int = 1): AR =
    {
        p = p_
        φ = durbinLevinson (p).slice (1, p+1)            // AR(p) coefficients: φ_0, ..., φ_p-1
        if (DEBUG) println (s"coefficients for AR($p): φ = $φ")
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector.
     */
    def parameters: VectoD = φ

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Apply the Durbin-Levinson Algorithm to iteratively compute the 'psi' matrix.
     *  The last row of the matrix gives 'AR' coefficients.
     *  @see www.stat.tamu.edu/~suhasini/teaching673/time_series.pdf
     */
    private def durbinLevinson: MatriD =
    {
        val psi = new MatrixD (ml+1, ml+1)
        val r   = new VectorD (ml+1); r(0) = ac(0)
        for (t <- 1 to ml) {
            var sum = 0.0
            for (j <- 1 until t) sum += psi(t-1, j) * ac(t-j)
            val a = (ac(t) - sum) / r(t-1)
            psi(t, t) = a
            for (j <- 1 until t) psi(t, j) = psi(t-1, j) - a * psi(t-1, t-j)
            r(t) = r(t-1) * (1.0 - a * a)
        } // for
        if (DEBUG) println ("psi = " + psi)
        pacf = psi.getDiag ().slice (1, ml+1)            // PACF is the diagonal of the psi matrix
        psi                                              // return the psi matrix
    } // durbinLevinson

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a vector that is the predictions of a 'p'th order Auto-Regressive 'AR(p)' model.
     */
    def predict (): VectoD =
    {
        val xp = new VectorD (n)                         // forecasts for all x
        for (t <- 0 until n) {
            var sum = 0.0
            for (j <- 0 until p if t-j > 0) sum += φ(j) * x(t-1-j)
            xp(t) = sum
        } // for
        e = x - xp
        xp + mu                                          // return the vector of predicted values
    } // predict

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Produce the multi-step forecast for AR models.
     *  @param steps  the number of steps to forecast, must be at least one.
     */
    def forecast (steps: Int = 1): VectoD =
    {
        val xf = x.slice (x.dim - p) ++ new VectorD (steps)

        for (t <- p until xf.dim) {                      // start at t = p (enough data and first value to forecast)
            var sum = 0.0
            for (j <- 0 until p) sum += φ(j) * xf(t-1-j)
            xf(t) = sum
        } // for
        xf.slice (p) + mu                                // return the vector of forecasts
    } // forecast

} // AR class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `AR` companion object provides factory methods for the `AR` class.
 */
object AR
{
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create an `AR` object from a matrix.
     *  @param ty  the matrix holding the time vector 't' and time series data 'y'
     */
    def apply (ty: MatriD): AR = new AR (ty.col(0), ty.col(1))

} // AR object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARTest` object is used to test the `AR` class.
 *  > runMain scalation.analytics.ARTest
 */
object ARTest extends App
{
    val ran = Random ()
    val n = 100
    val t = VectorD.range (0, n)
    val y = VectorD (for (i <- 0 until n) yield t(i) + 10.0 * ran.gen)

    val ts = new AR (t, y)              // time series data: t and y

    ts.plotFunc (ts.acf, "ACF")

    // Build AR(1) and AR(2) models for the time series data

    ts.train (1)
    println (s"φ_a = ${ts.parameters}")
    new Plot (t, y, ts.predict (), "Plot of y, ar(1) vs. t") //, true)

    ts.train (2)
    println (s"φ_b = ${ts.parameters}")
    new Plot (t, y, ts.predict (), "Plot of y, ar(2) vs. t") //, true)

    ts.plotFunc (ts.pacf, "PACF")

} // ARTest object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARTest2` object is used to test the `AR` class.
 *  > runMain scalation.analytics.ARTest2
 */
object ARTest2 extends App
{
    val noise = Normal (0.0, 1.0)
    val n = 20
    val t = VectorD.range (0, n)
    val y = VectorD (for (i <- 0 until n) yield i + noise.gen)

    println (s"y = $y")

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

    val ts = new AR (t, y)               // time series data: t and y

    // Build AR(1) model for the (differenced) time series data

    val φ_a = ts.train (p)
    println (s"φ_a = $φ_a")
    new Plot (t, y, ts.predict (), s"Plot of y, ar($p) vs. t") //, true)
    val ar_f = ts.forecast (steps)
    println (s"$steps-step ahead forecasts using AR($p) model = $ar_f")

} // ARTest2 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARTest3` object is used to test the `AR` class.
 *  Forecasting lake levels.
 *  @see cran.r-project.org/web/packages/fpp/fpp.pdf
 * > runMain scalation.analytics.ARTest3
 */
object ARTest3 extends App
{
    val t = VectorD.range (0, 98)
    val y = VectorD (580.38, 581.86, 580.97, 580.80, 579.79, 580.39, 580.42, 580.82, 581.40, 581.32,
                     581.44, 581.68, 581.17, 580.53, 580.01, 579.91, 579.14, 579.16, 579.55, 579.67,
                     578.44, 578.24, 579.10, 579.09, 579.35, 578.82, 579.32, 579.01, 579.00, 579.80,
                     579.83, 579.72, 579.89, 580.01, 579.37, 578.69, 578.19, 578.67, 579.55, 578.92,
                     578.09, 579.37, 580.13, 580.14, 579.51, 579.24, 578.66, 578.86, 578.05, 577.79,
                     576.75, 576.75, 577.82, 578.64, 580.58, 579.48, 577.38, 576.90, 576.94, 576.24,
                     576.84, 576.85, 576.90, 577.79, 578.18, 577.51, 577.23, 578.42, 579.61, 579.05,
                     579.26, 579.22, 579.38, 579.10, 577.95, 578.12, 579.75, 580.85, 580.41, 579.96,
                     579.61, 578.76, 578.18, 577.21, 577.13, 579.10, 578.25, 577.91, 576.89, 575.96,
                     576.80, 577.68, 578.38, 578.52, 579.74, 579.31, 579.89, 579.96)

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

    val ts = new AR (t, y)             // time series data: t and y

    // Build AR(2) model for the time series data

    val φ_a = ts.train (p)
    println (s"φ_a = $φ_a")
    new Plot (t, y, ts.predict (), s"Plot of y, ar($p) vs. t") //, true)
    val ar_f = ts.forecast (steps)
    println (s"$steps-step ahead forecasts using AR($p) model = $ar_f")

} // ARTest3 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `ARTest4` object is used to test the `AR` class.
 * > runMain scalation.analytics.ARTest4
 */
object ARTest4 extends App
{
    val path = BASE_DIR + "travelTime.csv"

    val data = MatrixD (path)

    val t = data.col(0)
    val y = data.col(1)

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

    val ts = new AR (t, y)               // time series data: t and y

    println (s"y = $y")

    // Build AR(1) model for the time series data

    val φ_a = ts.train (p)
    println (s"φ_a = $φ_a")
    new Plot (t, y, ts.predict (), s"Plot of y, ar($p) vs. t") //, true)
    val ar_f = ts.forecast (steps)
    println (s"$steps-step ahead forecasts using AR($p) model = $ar_f")

} // ARTest4 object