* 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. *------------------------------------------------------------------------------ * @param y_ the original input vector (time series data) * @param d the order of Integration/simple differencing * @param dd the order of seasonal differencing * @param period the seasonal period */ class SARIMA (y_ : VectoD, d: Int = 0, dd: Int = 0, period: Int = 1) extends ForecasterVec (y_, 1, null) // FIX degrees of freedom: dfm and df { private val DEBUG = true // debug flag private val differenced = d > 0 || dd > 0 // flag indicating whether differencing will be applied private var mu = -0.0 // sample mean (-0.0 means unassigned) private var μ = -0.0 // population mean estimated using MLE private var sig2 = -0.0 // sample variance private var σ2 = -0.0 // population variance estimated using MLE 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 params = 0 // number of parameters 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 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 δ: Double = -0.0 // drift/constant term init (y_) // initialize most of the above parameters if (DEBUG) { println ("z = " + z) // zero-mean time series println ("mu = " + mu) // mean println ("sig2 = " + sig2) // variance } // for //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 v the working vector/time-series */ def init (v: VectoD) { y = v 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 zp = new VectorD (z.dim) // vector of predicted/fitted values prior to undifferencing/un-zero mean e = new VectorD (z.dim) // vector of residuals sig2 = z.variance // the sample variance } // init //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return the model name including its current hyper-parameter. */ override def modelName: String = { if (period > 1) s"SARIMA ($p, $d, $q) x ($pp, $dd, $qq)_${period}" else s"SARIMA ($p, $d, $q)" } // modelName //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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): SARIMA = { p = p_ ; φ = new VectorD (p) // p and AR coefficients q = q_ ; θ = new VectorD (q) // q and MA coefficients pp = pp_; φφ = new VectorD (pp) // pp and Seasonal AR coefficients qq = qq_; θθ = new VectorD (qq) // qq and Seasonal MA coefficients cap = Array (p, q, pp*period, qq*period).max // greatest lag params = p + q + pp + qq + (if (differenced) 0 else 1) // number of parameters this } // setPQ //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Train/fit an `SARIMA` 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 (): SARIMA = { val solver = new QuasiNewton (nll) // nonlinear optimizer val b = new VectorD (params + 1) // parameter values 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 δ = μ * (1 - φ.sum) // update drift value // δ = stats.mu * (1 - φ.sum) if (DEBUG) { // showParameterEstimates () // show ARMA parameters println (s"θθ = $θθ") // seasonal moving average coeficients println (s"φφ = $φφ") // seasonal auto-regressive coeficients } // 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 input 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+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) 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 <- 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 } // for -ll (e.normSq / m, σ2, m) // return negative log likelihood } // updateFittedValues //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Return the vector of preducted/fitted values on the training/full data. */ override def predictAll (): VectoD = { if (differenced) transformBack (zp, z_, y, d, dd, period) else zp + μ } //predictAll override def predictAllz (): VectoD = predictAll () - μ //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** 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 } // parameter //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** Produce n-steps-ahead forecast for SARIMA 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 if tz-cap+i >= 0) { 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 = 0.0 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 - predictAll () 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 = predictAll ()): SARIMA = { resetDF (params, y.dim - params) 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): SARIMA = eval_ (y - yp, y, yp) */ // missing method def forwardSel (cols: Set [Int], index_q: Int): (Int, ForecasterVec) = ??? } // SARIMA class //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMATest` object is used to test the `SARIMA` class. * > runMain scalation.analytics.forecaster.SARIMATest */ object SARIMATest extends App { println ("SARIMA!") val m = 100 val noise = Uniform (-5, 5) val y = VectorD (for (i <- 0 until m) yield i + noise.gen) val ts = new SARIMA (y, 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 SARIMA(1, 0, 0), SARIMA(2, 0, 0) and SARIMA(0, 0, 3) models // ts.setPQ (1) // ts.train ().eval_ () // println (ts.report) // new Plot (null, y, ts.predictAll (), "Plot of y, SARIMA(1, 0, 0) vs. t", true) // // ts.setPQ (2) // ts.train ().eval_ () // println (ts.report) // new Plot (null, y, ts.predictAll (), "Plot of y, SARIMA(2, 0, 0) vs. t", true) ts.setPQ (0, 3) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.predictAll () , "Plot of y, SARIMA(0, 0, 3) vs. t", true) println (ts.forecastAll(3)) } // SARIMATest object //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMATest2` object is used to test the `SARIMA` class. * > runMain scalation.analytics.forecaster.SARIMATest2 */ object SARIMATest2 extends App { val m = 30 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 d = 1 // order of differencing val q = 1 val steps = 2 // number of steps for the forecasts val ts = new SARIMA (y, d) // time series data: y vs. t // Build ARIMA(p, 0), ARIMA(0, q) and ARIMA(p, q) models for the (differenced) time series data ts.setPQ (p) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.predictAll (), s"Plot of y, SARIMA($p, 0, 0) vs. t", true) val ar_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using SARIMA($p, 0, 0) model = $ar_f") ts.setPQ (0, q) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.predictAll (), s"Plot of y, SARIMA(0, 0, $q) vs. t", true) val ma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using SARIMA(0, 0, $q) model = $ma_f") ts.setPQ (p, q) ts.train ().eval_ () println (ts.report) new Plot (null, y, ts.predictAll (), s"Plot of y, SARIMA($p, 0, $q) vs. t", true) val arma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using SARIMA($p, 0, $q) model = $arma_f") } // SARIMATest2 object //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMATest3` object is used to test the `SARIMA` class. * Forecasting lake levels. * @see cran.r-project.org/web/packages/fpp/fpp.pdf * > runMain scalation.analytics.forecaster.SARIMATest3 */ object SARIMATest3 extends App { import ForecasterVec.y // val (pp, dd, qq) = (1, 1, 1) // seasonal hyper-parameters // val period = 4 // seasonal period // val d = 0 // apply d-th order differencing - no differencing val d = 1 // apply d-th order differencing - first differences for (h <- 1 to 1) { // forecasting horizon for (p <- 0 to 2) { // auto-regressive hyper-parameter settings for (q <- 0 to 2 if p+q > 0) { // moving-average hyper-parameter settings banner (s"Test3: SARIMA ($p, $d, $q) with h = $h") val mod = new SARIMA (y, d) // create an SARIMA model mod.setPQ (p, q).train ().eval_ () // set p and q, train and evaluate the SARIMA model println (mod.report) val yp = mod.predictAll () // predicted values val yf = mod.forecastAll (h)(h) // forecasted values - h steps ahead // val yf2 = mod.forecastAll2 (h).col(h) // forecasted values - h steps ahead new Plot (null, y, yp, s"Plot of y & yp, predicted SARIMA ($p, $d, $q) vs. t", true) new Plot (null, y, yf, s"Plot of y & yf, forecasted (h = $h) SARIMA ($p, $d, $q) vs. t", true) // new Plot (null, y, yf2, s"Plot of y & yf2, forecasted2 (h = $h) SARIMA ($p, $d, $q) vs. t", true) if (h == 1) ForecasterVec.differ (yp, yf, allow = true) } // for } // for } // for } // SARIMATest3 object //:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: /** The `SARIMATest4` object is used to test the `SARIMA` class. * > runMain scalation.analytics.forecaster.SARIMATest4 */ object SARIMATest4 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 SARIMA (y, d) // time series data: y vs. t println (s"y = $y") // Build SARIMA(p, d, q) models ts.setPQ (p) ts.train ().eval_ () println (ts.report) new Plot (t, y, ts.predictAll (), s"Plot of y, SARIMA($p, $d, 0) vs. t", true) val ar_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using SARIMA($p, $d, $q) model = $ar_f") ts.setPQ (0, q) ts.train ().eval_ () println (ts.report) new Plot (t, y, ts.predictAll (), s"Plot of y, SARIMA(0, $d, $q) vs. t", true) val ma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using SARIMA($p, $d, $q) model = $ma_f") ts.setPQ (p, q) ts.train ().eval_ () println (ts.report) new Plot (t, y, ts.predictAll (), s"Plot of y, SARIMA($p, $d, $q) vs. t", true) val arma_f = ts.forecast (h = steps) println (s"$steps-step ahead forecasts using SARIMA($p, $d, $q) model = $arma_f") } // SARIMATest4 object