//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John Miller
 *  @version 1.6
 *  @date    Sat Feb 15 15:44:16 EST 2020
 *  @see     LICENSE (MIT style license file).
 *
 *  @title   Model: Moving-Average
 *
 *  @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     https://maryclare.github.io/atsa/content/notes/notes_3.pdf
 *  @see     https://math.unice.fr/~frapetti/CorsoP/Chapitre_4_IMEA_1.pdf
 */

package scalation.analytics
package forecaster

import scala.collection.mutable.Set
import scala.math.{abs, sqrt}

import scalation.linalgebra.{MatriD, MatrixD, VectoD, VectorD}
import scalation.math.{double_exp, noDouble}
import scalation.minima.GoldenSectionLS
import scalation.plot.Plot
import scalation.random.{Normal, Random, Uniform}
import scalation.stat.{Statistic, StatVector, vectorD2StatVector}
import scalation.util.banner

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

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MA_1` class provides basic time series analysis capabilities for 'MA_1'
 *  models.  In an 'MA_1' model, 1 refers to the order of the Moving-Average
 *  component of the model.  `MA_1` 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 the noise/residuals/shocks:
 *  <p>
 *      y_t = c + θ e_t-i + e_t
 *  <p>
 *  where 'c' is a constant, 'θ' is the moving-average parameters vector,
 *  and 'e' is the noise vector.
 *------------------------------------------------------------------------------
 *  @param y       the response vector (time series data)
 *  @param hparam  the hyper-parameters
 */
class MA_1 (y: VectoD, hparam: HyperParameter = null)
      extends ForecasterVec (y, 1, hparam) with NoFeatureSelectionF
{
    private val DEBUG = true                                    // debug flag

    private var θ:  VectoD = null                               // MA(1) parameter
    e = new VectorD (y.dim)                                     // residual/error vector

//  private val z    = y - mu                                   // work with mean zero time series
    private var zp: VectoD = null                               // centered predicted values
//  private val (sig2, acv, acr) = acf (z, ml)                  // variance, auto-covariance, auto-correlation

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Train an `MA_1` model on the time-series data in vector 'y_'.
     *  Estimate the parameters vector 'θ' for a 'q'th order a Moving-Average 'MA(1)'
     *  model.
     *  <p>
     *      z_t = θ_0 * e_t-1 + e_t
     *  <p>
     *  @param x_null  the data/input vector (ignored)
     *  @param y_      the response/output vector (currently only works for y)
     */
    override def train (x_null: MatriD, y_ : VectoD): MA_1 =
    {
        super.train (null, y_)
        val ρ = stats.acr(1)
//      θ = VectorD (directEstimate (ρ))                       // MoM
        θ = VectorD (iterativeEstimate)                        // LSE
        if (DEBUG) println (s"train: ρ = $ρ, θ = $θ")
        this
    } // train

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector (θ) by inverting the autocorrelation function.
     *  <p>
     *      -θ / (1 + θ^2) = ρ
     *  <p>
     *  It uses the Method of Moments (less accurate).
     *  @see equation 3.3.7 in Box-Jenkins
     *  @see people.stat.sc.edu/hitchcock/stat520ch7slides.pdf
     *  @param ρ  the first order auto-correlation
     */
    private def directEstimate (ρ: Double): Double =
    {
        if (abs (ρ) > 0.5) flaw ("estimate", s"the first order auto-correlation $ρ is outside the invertible interval")
        val r2 = 2 * ρ
        val rt = sqrt (1 - 4 * ρ~^2)
        val θ = (-1 + rt) / r2
        if (abs (θ) < 1) θ else (-1 - rt) / r2
    } // directEstimate

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Given a value for parameter 'θ', compute the Conditional Sum of Squared Errors.
     *  @see people.stat.sc.edu/hitchcock/stat520ch7slides.pdf
     *  @param θ  the given θ parameters
     */
    def csse (θ: Double): Double =
    {
        e(0)      = 0.0
        var sumSq = 0.0
        for (t <- 1 until e.dim) {
            e(t) = z(t) - θ * e(t-1)
            sumSq += e(t) ~^ 2
        } // for
        sumSq
    } // csse

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector (θ) by minimizing the Conditional Sum of Squared Errors.
     *  Least Squares Estimation.
     */
    private def iterativeEstimate: Double =
    {
        val golds = new GoldenSectionLS (csse _)
        val θ = golds.lsearch (1.0, -1.0)                               // complete range
        println (s"golds: θ = $θ, csse (θ) = ${csse (θ)}")
        θ
    } // iterativeEstimate

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a detailed summary of the trained model.
     */
    def summary: String =
    {
        super.summary ("θ", s"${stats.mu} + $θ dot [e_(t-1)]")
    } // summary

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return the parameter vector (θ).
     */
    def parameter: VectoD = θ

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Return a vector of centered predictions of an MA_1 model and update the residuals.
     */
    def predictAllz (): VectoD =
    {
        val zp = new VectorD (m)
        zp(0)  = z(0)
        for (t <- 1 until m) {
            zp(t) = θ(0) * e(t-1)
            e(t)  = z(t) - zp(t)
        } // for
        zp
    } // predictAllz

    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** 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).
     *  @param h  the maximum forecasting horizon, number of steps ahead to produce forecasts
     */
    def forecastAll (h: Int): MatriD =
    {
        yf     = new MatrixD (m, h+1)                                    // forecasts for all time points t & horizons to h
        val ef = new VectorD (m)                                         // error vector based on forecasts
        yf.setCol (0, y)                                                 // first column is actual values, horizon 0
        for (k <- 1 to h) {
            yf(0, k) = y(0)                                              // copy first actual value
            ef(0)    = 0.0                                               // copy => no error at first
            for (t <- 1 until m) {                                       // forecast the rest
                 yf(t, k) = stats.mu + θ(0) * ef(t-1)
                 ef(t)    = y(t) - yf(t, k)
            } // for
            if (DEBUG) println (s"forecastAll: yf.col ($k) = ${yf.col (k)}")
        } // for
        yf                                                               // return matrix of forecasted values
    } // forecastAll
 
} // MA_1 class


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MA_1` companion object provides factory methods for the `MA_1` class.
 */
object MA_1
{
    //::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
    /** Create an `MA_1` object.
     *  @param y  the response vector (time series data)
     */
    def apply (y: VectoD): MA_1 =
    {
        new MA_1 (y, null)
    } // apply

} // MA_1 object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MA_1Test` object is used to test the `MA_1` class.
 *  <p>
 *      y_t = c + θ e_{t-1} + e_t
 *  <p>
 *  > runMain scalation.analytics.forecaster.MA_1Test
 */
object MA_1Test extends App
{
    val m     = 100
    val noise = Normal (0, 100)
//  val noise = Uniform (-30, 30)
    val y = new VectorD (m)
    var e = noise.gen
    for (i <- y.range) {
        val ee = noise.gen
        y(i) = 0.5 * e + ee 
        e = ee
    } // for

    val ma = new MA_1 (y)                                  // time series model MA_1
    val ar = new AR (y)                                    // time series model AR
    
    banner ("Build AR(1) model")
    ar.train (null, y).eval ()
    val yp = ar.predictAll ()
    println (ar.report)
    new Plot (null, y, yp, "Plot of y, AR(1) vs. t", true)

    banner ("Build AR(2) model")
    ar.retrain (2).eval ()
    val yp2 = ar.predictAll ()
    println (ar.report)
    new Plot (null, y, yp2, "Plot of y, AR(2) vs. t", true)

    banner ("Build MA_1 model")
    ma.train (null, y).eval ()
    val yp3 = ma.predictAll ()
    println (ma.report)
    new Plot (null, y, yp3, "Plot of y, MA(1) vs. t", true)

    ar.plotFunc2 (ar.acF, "ACF")
    ar.plotFunc2 (ar.pacF, "PACF")

} // MA_1Test object


//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** The `MA_1Test2` object is used to test the `MA_1` class on simulated data with
 *  normally distributed errors.
 *  > runMain scalation.analytics.forecaster.MA_1Test2
 */
object MA_1Test2 extends App
{
    println ("TBD")

} // MA_1Test2 object


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

    var ma: MA_1 = null
    for (h <- 1 to 2) {                                            // forecasting horizon
        banner (s"Test3: MA_1 with h = $h")
        ma = new MA_1 (y)                                          // MA_1 model for time series data
        ma.train (null, y).eval ()                                 // train the model and evaluate
        val yf  = ma.forecastAll (h)
        println (s"yf = $yf")
        val yf2 = testInSamp (y, 1, ma, h)                         // in-sample test
        assert (yf == yf2)
        println (ma.summary)
    } // for

    banner ("Select model based on ACF and PACF")
    ma.plotFunc2 (ma.acF, "ACF")
    ma.plotFunc2 (ma.pacF, "PACF")

    banner ("Rolling Validation")
    val stats = SimpleRollingValidation.crossValidate2 (ma, kt_ = 5)
//  val stats = RollingValidation.crossValidate2 (ma, kt_ = 5)
    Fit.showQofStatTable (stats)

} // MA_1Test3 object