Merge pull request #41 from josepablocam/ljungbox

Ljung-Box test

Author: sryza

src/main/scala/com/cloudera/sparkts/TimeSeriesStatisticalTests.scala

@@ -272,7 +272,27 @@ object TimeSeriesStatisticalTests {
     val bgstat = nObs * auxOLS.calculateRSquared()
     (bgstat, 1 - new ChiSquaredDistribution(maxLag).cumulativeProbability(bgstat))
   }
-  
+
+  /**
+   * Ljung-Box test for serial correlation in residuals up to lag `maxLag`. The null hypothesis
+   * is that values are independently distributed up to the given lag. The alternate hypothesis
+   * is that serial correlation is present. The test statistic follows a Chi-Squared distribution
+   * with `maxLag` degrees of freedom. See [[https://en.wikipedia.org/wiki/Ljung%E2%80%93Box_test]]
+   * for more information.
+   * @param residuals
+   * @return the test statistic and the p-value associated with it.
+   */
+  def lbtest(residuals: Vector[Double], maxLag: Int): (Double, Double) = {
+    val autoCorrs = UnivariateTimeSeries.autocorr(residuals, maxLag)
+    val n = residuals.length
+    val adjAutoCorrs = autoCorrs.toArray.zipWithIndex.map { case (p, k) =>
+      (p * p) / (n - k - 1)
+    }
+    val testStatistic = n * (n + 2) * adjAutoCorrs.sum
+    val pValue = 1 - new ChiSquaredDistribution(maxLag).cumulativeProbability(testStatistic)
+    (testStatistic, pValue)
+  }
+
   /**
    * Breusch-Pagan test for heteroskedasticity in a model
    * The statistic follows a X^2 distribution with (# of regressors - 1) degrees of freedom

src/test/scala/com/cloudera/sparkts/TimeSeriesStatisticalTestsSuite.scala

@@ -86,4 +86,20 @@ class TimeSeriesStatisticalTestsSuite extends FunSuite with ShouldMatchers {
     // there should be evidence of heteroskedasticity
     bptest(new DenseVector(resids2), new DenseMatrix(x.length, 1, x))._2 should be < pthreshold
   }
+
+  test("ljung-box test") {
+    val rand = new MersenneTwister(5L)
+    val n = 100
+    val indep = Array.fill(n)(rand.nextGaussian)
+    val vecIndep = new DenseVector(indep)
+    val (stat1, pval1) = lbtest(vecIndep, 1)
+    pval1 > 0.05
+
+    // serially correlated
+    val coef = 0.3
+    val dep = indep.scanLeft(0.0) { case (prior, curr) => prior * coef + curr }.tail
+    val vecDep = new DenseVector(dep)
+    val (stat2, pval2) = lbtest(vecDep, 2)
+    pval2 < 0.05
+  }
 }