* y_t = c + Σ(φ_i y_t-i) + Σ(θ_i e_t-i) + βχ + e_t *
* where 'c' is a constant, 'φ' is the autoregressive coefficient vector, * 'θ' is the moving-average coefficient vector, and 'e' is the noise vector. * If 'd' > 0, then the time series must be differenced first before applying * the above model. Seasonal differencing, autoregressive and moving average * factors can be incorporated into the model by applying seasonal differencing * (possibly in addition to simple differencing) first, then add the seasonal * autoregressive and moving average terms, that rely on lagged values and errors, * respectively, from one or more seasonal periods in the past, on the right hand * side of the equation. Exogenous/External regressor may also be added to the * right-hand size of the model in a similar manner to Regression models. *------------------------------------------------------------------------------ * @param y_ the original input vector (time series data) * @param χ the exogenous vector (time series data) * @param lg the lag for the exogenous variable, note need lg >= horizon h * @param d the order of Integration/simple differencing * @param dd the order of seasonal differencing * @param period the seasonal period */ class SARIMAX (y_ : VectoD, χ: VectoD, lg: Int, d: Int = 0, dd: Int = 0, period: Int = 1) extends ForecasterVec (y_, 1, null) // FIX degrees of freedom: dfm and df { private val DEBUG = false // debug flag private val differenced = d > 0 || dd > 0 // whether differencing will be applied private var m_ = 0 // size of the time series (after differencing) private var mu = 0.0 // the sample mean private var μ = 0.0 // the true population mean (estimated using MLE) private var sig2 = 0.0 // the sample variance private var σ2 = 0.0 // the true population variance (estimated using MLE) private var y: VectoD = null // current y, initially set to y_ private var z_ : VectoD = null // intermediate results after simple differencing private var zp: VectoD = null // vector of predicted/fitted values prior to undifferencing/un-zero mean setTS (y_) // initialize most of the above parameters private var p = 0 // AR order private var q = 0 // MA order private var pp = 0 // seasonal AR order private var qq = 0 // seasonal MA order private var cap = 0 // max among p, q, pp*period, and qq*period private var φ: VectoD = null // AR(p) coefficients private var φφ: VectoD = null // seasonal AR(pp) coefficients private var θ: VectoD = null // MA(q) coefficients private var θθ: VectoD = null // seasonal MA(qq) coefficients private var β = 0.0 // exogeneous variable coefficient if (DEBUG) { println ("m_ = " + m_) // number of data points println ("z = " + z) // zero-mean time series println ("mu = " + mu) // mean println ("sig2 = " + sig2) // variance } // for //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return the model name including its current hyper-parameter. */ override def modelName: String = { if (period > 1) s"SARIMAX ($p, $d, $q) x ($pp, $dd, $qq)_${period}" else s"SARIMAX ($p, $d, $q)" } // modelName //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set/change the internal time series. May be used to set the time series * to a different time window (typically future when new data become available) * in order to produce newer forecast (typically with the new data) without re-training * the model for parameters (use existing parameters from previous training). * @param ts the new time series */ def setTS (ts: VectoD): Unit = { y = ts mu = y.mean // the sample mean val (z0, z0_) = difference (y, d, dd, period) // difference the time series z = if (differenced) z0 else z0.copy // time series after differencing // deep copy needed for undifferenced time series to avoid modifying original y z_ = z0_ // intermediate results after simple differencing m_ = z.dim // size of the time series (after differencing) e = new VectorD (m_) // vector of residuals zp = new VectorD (m_) // vector of predicted/fitted values prior to undifferencing/un-zero mean sig2 = z.variance // the sample variance } // setTS //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Set values for 'p', 'q', 'pp' and 'qq'. * @param p_ the order of the AR part of the model * @param q_ the order of the MA part of the model * @param pp_ the order of the Seasonal AR part of the model * @param qq_ the order of the Seasonal MA part of the model */ def setPQ (p_ : Int = 0, q_ : Int = 0, pp_ : Int = 0, qq_ : Int = 0): Unit = { p = p_ ; if (p > 0) φ = new VectorD (p) q = q_ ; if (q > 0) θ = new VectorD (q) pp = pp_; if (pp > 0) φφ = new VectorD (pp) qq = qq_; if (qq > 0) θθ = new VectorD (qq) cap = Array (p, q, pp*period, qq*period).max } // setPQ //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Train/fit an `SARIMAX` model to the times series data. Must call 'SetPQ' first. */ def train (): SARIMAX = { val solver = new QuasiNewton (nll) // nonlinear optimizer val b = new VectorD (p+pp+q+qq+2+(if (differenced) 0 else 1)) // parameter values b(b.size-3) = β // exogeneous variable coefficient if (! differenced) b(b.size-2) = mu // sample mean, initial est. for μ parameter b(b.size-1) = sqrt (sig2) // sample standard deviation, initial est. for σ parameter solver.solve (b) // find b the maximizes likelihood if (DEBUG) { println (s"φ = $φ") // autoregressive coefficients println (s"φφ = $φφ") // seasonal autoregressive coefficients println (s"θ = $θ") // moving average coefficients println (s"θθ = $θθ") // seasonal moving averag coefficients println (s"β = $β") // exogeneous variable coefficient println (s"mu = $mu") // sample mean println (s"μ = $μ") // mean parameter in MLE println (s"σ2 = $σ2") // variance parameter in MLE } // if this } // train //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The negative log-likelihood function to be minimized. * @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 != p+pp+q+qq+2+(if (differenced) 0 else 1)) flaw ("nll: ", "input param vector size incorrect") for (i <- 0 until p) φ(i) = b(i) for (i <- p until p+pp) φφ(i-p) = b(i) for (i <- p+pp until p+pp+q) θ(i-p-pp) = b(i) for (i <- p+pp+q until p+pp+q+qq) θθ(i-p-pp-q) = b(i) β = b(b.size-3) if (! differenced) μ = b(b.size-2) σ2 = b.last~^2 updateFittedValues () } // nll //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Update 'xp', the vector of fitted values; 'e', the vector of errors; * ll, aic, aicc and bic. */ def updateFittedValues (): Double = { if (!differenced) for (i <- z.range) z(i) = y(i) - μ // for undifferenced time series, perform zero-mean using currently estimated μ for (t <- 0 until zp.dim) { var sum = β * χ(max (0, t - lg)) for (j <- 0 until p if t-j > 0) sum += φ(j) * z(t-1-j) for (j <- 0 until pp if t-(1+j)*period >= 0) sum += φφ(j) * z(t-(1+j)*period) for (j <- 0 until q if t-j > 0) sum += θ(j) * e(t-1-j) for (j <- 0 until qq if t-(1+j)*period >= 0) sum += θθ(j) * e(t-(1+j)*period) zp(t) = sum e(t) = z(t) - sum } // for -ll (e.normSq / m, σ2, m) } // updateFittedValues //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return the vector of fitted values on the training data. */ def fittedValues (): VectoD = if (differenced) transformBack (zp, z_, y, d, dd, period) else zp + μ override def predictAllz (): VectoD = fittedValues () - μ override def predictAll (): VectoD = fittedValues () //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return the parameter vector. */ override def parameter: VectoD = { var b = new VectorD (0) if (p > 0) b = b ++ φ if (q > 0) b = b ++ θ if (pp > 0) b = b ++ φφ if (qq > 0) b = b ++ θθ b ++ VectorD (β) } // parameter //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Produce n-steps-ahead forecast for SARIMAX 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 tz = t - d - dd * period // scale t to match vector z and e if (tz < cap) flaw ("forecast", s"tz ($tz) must be at least cap ($cap)") val zf = new VectorD (cap + h) val ef = new VectorD (cap + h) for (i <- 0 until cap) { zf(i) = z(tz-cap+i) ef(i) = e(tz-cap+i) } // for for (i <- cap until zf.dim) { // start at t = cap (enough data and first value to forecast) var sum = β * χ(max (0, t - lg)) for (j <- 0 until p) sum += φ(j) * zf(i-1-j) for (j <- 0 until pp) sum += φφ(j) * zf(i-(1+j)*period) for (j <- 0 until q) sum += θ(j) * ef(i-1-j) for (j <- 0 until qq) sum += θθ(j) * ef(i-(1+j)*period) zf(i) = sum } // for val f = zf.slice (cap) if (differenced) transformBackF (f, z_, y, d, dd, period, t) else 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 */ def forecastAll (h: Int = 1): MatriD = { val yf = new MatrixD (h+1, y.dim) // forecasts for all horizons h & time points t yf(0) = y // first row is actual values val st = cap + d + dd * period // define the starting index for (t <- st until y.dim) { 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(k, t+k-1) = ft(k-1) // place the forecasted values diagonally down into yf // (to align with the appropriate forecasting horizon/step) } // for } // for // fill in the blank values in the first few columns where no forecasts can be produced by copying values from previous rows for (k <- 1 to h; t <- 0 until st) yf(k, t) = y(t) // copy observed values from y for (k <- 2 to h; t <- st until st+k-1) yf(k, t) = yf(k-1, t) // copy forecasted values if (DEBUG) for (k <- 1 to h) println (s"for k = $k, sse = ${(y - yf(k)).normSq}") yf // return matrix of forecasted values } // forecastAll //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Obtain residuals in the original scale */ def residuals () = if (differenced) y - fittedValues () else 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 = fittedValues ()): SARIMAX = { val numParam = p+pp+q+qq+(if (differenced) 0 else 1) // intercept/mean counts as a parameter for undifferenced time series resetDF (numParam, y.dim-numParam) diagnose (e, y, yp) this } // eval_ //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Compute the error (difference between actual and predicted) and useful * diagnostics for the dataset. * @param y vector of observed values * @param yp vector of predicted values */ def eval_ (y: VectoD, yp: VectoD): SARIMAX = eval_ (y - yp, y, yp) // missing method def forwardSel (cols: Set [Int], index_q: Int): (Int, ForecasterVec) = ??? } // SARIMAX class //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMAXTest` object is used to test the `SARIMAX` class. * > runMain scalation.analytics.forecaster.SARIMAXTest */ object SARIMAXTest extends App { println ("SARIMAX!") val m = 100 val noise = Uniform (-5, 5) val y = VectorD (for (i <- 0 until m) yield i + noise.gen) val x = VectorD (for (i <- 0 until m) yield i / 2.0) val ts = new SARIMAX (y, x, 1, 1) // time series data: y vs. t, apply 1st order differencing // ts.plotFunc2 (ts.acf, "ACF") // must turn on DEBUG so that ACF is actually computed // Build ARIMAX(p, d, q) models // ts.setPQ (1) // ts.train () // new Plot (t, y, ts.fittedValues (), "Plot of y, ARIMAX(1, 0, 0) vs. t", true) // // ts.setPQ (2) // ts.train () // new Plot (t, y, ts.fittedValues (), "Plot of y, ARIMAX(2, 0, 0) vs. t", true) ts.setPQ (0, 3) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.fittedValues () , "Plot of y, ARIMAX(0, 0, 3) vs. t", true) println (ts.forecastAll(3)) } // SARIMAXTest object //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMAXTest2` object is used to test the `SARIMAX` class. * > runMain scalation.analytics.forecaster.SARIMAXTest2 */ object SARIMAXTest2 extends App { val m = 30 val noise = Normal (0, 2) val y = VectorD (for (i <- 0 until m) yield i + noise.gen) val x = VectorD (for (i <- 0 until m) yield i / 2.0) println (s"y = $y") val p = 1 val d = 1 // apply 1st order differencing val q = 1 val steps = 2 // number of steps for the forecasts val ts = new SARIMAX (y, x, 1, d) // time series data: y vs. t // Build ARIMAX(p, d, q) models ts.setPQ (p) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.fittedValues (), s"Plot of y, ARIMAX($p, $d, 0) vs. t", true) val ar_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX($p, $d, 0) model = $ar_f") ts.setPQ (0, q) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.fittedValues (), s"Plot of y, ARIMAX(0, $d, $q) vs. t", true) val ma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX(0, $d, $q) model = $ma_f") ts.setPQ (p, q) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.fittedValues (), s"Plot of y, ARIMAX($p, $d, $q) vs. t", true) val arma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX($p, $d, $q) model = $arma_f") } // SARIMAXTest2 object //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMAXTest3` object is used to test the `SARIMAX` class. * Forecasting lake levels. * @see cran.r-project.org/web/packages/fpp/fpp.pdf * > runMain scalation.analytics.forecaster.SARIMAXTest3 */ object SARIMAXTest3 extends App { import ForecasterVec.y val x = VectorD (for (i <- y.range) yield (y(max (0, i-1)) + y(i)) / 2.0) val p = 4 val q = 2 val d = 1 // order of differencing // val pp = 1 // val qq = 1 // val dd = 1 // val period = 4 val steps = 3 // number of steps for the forecasts val ts = new SARIMAX (y, x, 1, d) //, dd, period) // time series data: y vs. t // Build ARIMAX(p, d, q) models ts.setPQ (p, q) //, pp, qq) ts.train ().eval_ () println (ts.report) val xp = ts.fittedValues () println (s"xp = $xp") new Plot (null, y, xp, s"Plot of y, ARIMAX($p, $d, $q) vs. t", true) val arima_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX($p, $d, $q) model = $arima_f") println (ts.forecastAll(3)) } // SARIMAXTest3 object //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMAXTest4` object is used to test the `SARIMAX` class. * > runMain scalation.analytics.forecaster.SARIMAXTest4 */ object SARIMAXTest4 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 q = 1 val d = 1 // apply 1st order differencing val steps = 1 // number of steps for the forecasts val ts = new SARIMAX (y, t, 1, d) // time series data: y vs. t println (s"y = $y") // Build ARIMAX(p, d, q) models ts.setPQ (p) ts.train ().eval_ () println (ts.report) new Plot (t, y, ts.fittedValues (), s"Plot of y, ARIMAX($p, $d, 0) vs. t", true) val ar_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX($p, $d, 0) model = $ar_f") ts.setPQ (0, q) ts.train ().eval_ () println (ts.report) new Plot (t, y, ts.fittedValues (), s"Plot of y, ARIMAX(0, $d, $q) vs. t", true) val ma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX(0, $d, $q) model = $ma_f") ts.setPQ (p, q) ts.train ().eval_ () println (ts.report) new Plot (t, y, ts.fittedValues (), s"Plot of y, ARIMAX($p, $d, $q) vs. t", true) val arma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using ARIMAX($p, $d, $q) model = $arma_f") } // SARIMAXTest4 object