Merge pull request #416 from palatej/develop

Correction in CanovaHansen (seas and td). Code cleanup

jdplus-main-base/jdplus-sa-base-parent/jdplus-sa-base-core/src/main/java/jdplus/sa/base/core/tests/CanovaHansen.java

@@ -28,7 +28,12 @@
 import jdplus.toolkit.base.core.stats.RobustCovarianceComputer;
 import jdplus.toolkit.base.core.modelling.regression.PeriodicDummiesFactory;
 import jdplus.toolkit.base.api.data.DoubleSeq;
+import jdplus.toolkit.base.api.stats.StatisticalTest;
+import jdplus.toolkit.base.api.stats.TestType;
+import jdplus.toolkit.base.core.data.DataBlock;
+import jdplus.toolkit.base.core.dstats.F;
 import jdplus.toolkit.base.core.math.matrices.FastMatrix;
+import jdplus.toolkit.base.core.stats.tests.TestsUtility;
 
 /**
  *
@@ -103,7 +108,7 @@ public Builder startPosition(int startPosition) {
  public CanovaHansen build() {
  FastMatrix x = sx();
  LinearModel lm = buildModel(x);
- return new CanovaHansen(x, lm, winFunction, truncationLag);
+ return new CanovaHansen(lm, x.getColumnsCount(), winFunction, truncationLag);
  }
 
  private FastMatrix sx() {
@@ -114,17 +119,18 @@ private FastMatrix sx() {
  --len;
  }
  switch (type) {
- case Dummy: {
+ case Dummy -> {
  PeriodicDummies vars = new PeriodicDummies((int) period);
  return PeriodicDummiesFactory.matrix(vars, len, pos);
  }
- case Trigonometric: {
+ case Trigonometric -> {
  TrigonometricSeries vars = TrigonometricSeries.regular((int) period);
  return vars.matrix(len, pos);
  }
- default:
+ default -> {
  TrigonometricSeries vars = TrigonometricSeries.all(period, nh);
  return vars.matrix(len, pos);
+ }
  }
 
  }
@@ -139,14 +145,15 @@ private LinearModel buildModel(FastMatrix sx) {
  builder.y(s);
  }
  switch (type) {
- case Dummy ->  {
+ case Dummy -> {
  builder.addX(sx);
  }
- case Trigonometric ->  {
+ case Trigonometric -> {
  builder.addX(sx)
  .meanCorrection(true);
  }
- default -> builder.addX(sx)
+ default ->
+ builder.addX(sx)
  .meanCorrection(true);
  }
  return builder.build();
@@ -160,39 +167,57 @@ public DoubleSeq getE() {
  return u;
  }
 
- private final FastMatrix x, xe, cxe, omega;
+ private final FastMatrix x, xe, cxe, phi;
  private final DoubleSeq c, u;
+ private final int nx, ns;
 
- private CanovaHansen(final FastMatrix x, final LinearModel lm, final WindowFunction winFunction, int truncationLag) {
- this.x = x;
+ private CanovaHansen(final LinearModel lm, int ns, final WindowFunction winFunction, int truncationLag) {
+ this.ns=ns;
+ x = lm.variables();
+ nx=x.getColumnsCount();
  LeastSquaresResults olsResults = Ols.compute(lm);
- c=olsResults.getCoefficients();
+ c = olsResults.getCoefficients();
  u = lm.calcResiduals(c);
- xe = x.deepClone();
  // multiply the columns of x by e
+ xe=x.deepClone();
  xe.applyByColumns(col -> col.apply(u, (a, b) -> a * b));
- omega = RobustCovarianceComputer.covariance(xe, winFunction, truncationLag);
+ phi = RobustCovarianceComputer.covariance(xe, winFunction, truncationLag);
  cxe = xe.deepClone();
  cxe.applyByColumns(col -> col.cumul());
  }
 
  public double test(int var) {
- return computeStat(omega.extract(var, 1, var, 1), cxe.extract(0, cxe.getRowsCount(), var, 1));
+ int dx=nx-ns, dvar=var+dx;
+ return computeStat(phi.extract(dvar, 1, dvar, 1), cxe.extract(0, cxe.getRowsCount(), dvar, 1));
  }
 
  public double test(int var, int nvars) {
- return computeStat(omega.extract(var, nvars, var, nvars), cxe.extract(0, cxe.getRowsCount(), var, nvars));
+ int dx=nx-ns, dvar=var+dx;
+ return computeStat(phi.extract(dvar, nvars, dvar, nvars), cxe.extract(0, cxe.getRowsCount(), dvar, nvars));
  }
 
  public double testAll() {
- return computeStat(omega, cxe);
+ return test(0, ns);
+ }
+
+ public StatisticalTest seasonalityTest() {
+ int dx=nx-ns;
+ FastMatrix rcov = robustCovarianceOfCoefficients().extract(dx, ns, dx, ns);
+ SymmetricMatrix.lcholesky(rcov);
+ LowerTriangularMatrix.toLower(rcov);
+ DataBlock b = DataBlock.of(c.extract(dx, ns));
+ LowerTriangularMatrix.solveLx(rcov, b);
+ double fval = b.ssq() / ns;
+ F f = new F(ns, x.getRowsCount() - c.length());
+ return TestsUtility.testOf(fval, f, TestType.Upper);
  }
 
  private double computeStat(FastMatrix O, FastMatrix cx) {
- int n = cx.getRowsCount(), nx = cx.getColumnsCount();
+ int n = cx.getRowsCount(), ncx = cx.getColumnsCount();
  // compute tr( O^-1*xe'*xe)
  // cusum
- FastMatrix FF = FastMatrix.square(nx);
+ // FF = X'X
+ FastMatrix FF = FastMatrix.square(ncx);
  for (int i = 0; i < n; ++i) {
  FF.addXaXt(1, cx.row(i));
  }
@@ -206,18 +231,19 @@ private double computeStat(FastMatrix O, FastMatrix cx) {
  return tr / (n * n);
  }
 
-// private FastMatrix robustCovarianceOfCoefficients() {
-// FastMatrix Lo = omega.deepClone();
-// SymmetricMatrix.lcholesky(Lo);
-//
-// FastMatrix Lx = SymmetricMatrix.XtX(x);
-// SymmetricMatrix.lcholesky(Lx);
-// LowerTriangularMatrix.solveLX(Lx, Lo);
-// LowerTriangularMatrix.solveLtX(Lx, Lo);
-//
-// FastMatrix XXt = SymmetricMatrix.XXt(Lo);
-// XXt.mul(xe.getRowsCount());
-// return XXt;
-// }
+ private FastMatrix robustCovarianceOfCoefficients() {
+ FastMatrix Lo = phi.deepClone();
+ SymmetricMatrix.lcholesky(Lo);
+ LowerTriangularMatrix.toLower(Lo);
+
+ FastMatrix Lx = SymmetricMatrix.XtX(x);
+ SymmetricMatrix.lcholesky(Lx);
+ LowerTriangularMatrix.solveLX(Lx, Lo);
+ LowerTriangularMatrix.solveLtX(Lx, Lo);
+
+ FastMatrix XXt = SymmetricMatrix.XXt(Lo);
+ XXt.mul(xe.getRowsCount());
+ return XXt;
+ }
 
 }

jdplus-main-base/jdplus-sa-base-parent/jdplus-sa-base-core/src/test/java/jdplus/sa/base/core/tests/CanovaHansenTest.java

@@ -39,6 +39,7 @@ public void testUnempl_dummy() {
  .dummies(4)
  .truncationLag(4)
  .build();
+// System.out.println(ch.seasonalityTest());
 
 // for (int i = 0; i < 4; ++i) {
 // System.out.println(ch.test(i));
@@ -53,23 +54,24 @@ public void testUnempl_trig() {
  .trigonometric(4)
  .truncationLag(4)
  .build();
+// System.out.println(ch.seasonalityTest());
 // System.out.println(ch.test(0, 2));
 // System.out.println(ch.test(2));
 // System.out.println(ch.testAll());
  }
 
  @Test
  public void testP_dummy() {
- DoubleSeq y = delta(log(Data.TS_PROD), 1).getValues();
+ DoubleSeq y = log(Data.TS_PROD).getValues();
 // System.out.println("dummies");
  CanovaHansen ch = CanovaHansen.test(y)
  .dummies(12)
- .lag1(false)
- .truncationLag(12)
+ .lag1(true)
+ .truncationLag(15)
  .startPosition(1)
  .build();
- //  System.out.println(ch.robustTestCoefficients());
-//
+// System.out.println(ch.seasonalityTest());
+
 // for (int i = 0; i < 12; ++i) {
 // System.out.println(ch.test(i));
 // }
@@ -93,7 +95,7 @@ public void testP_trig() {
 // }
  // System.out.println(ch.robustTestCoefficients());
 // System.out.println(all);
- }
+  }
 
  @Test
  public void testW_trig() {
@@ -108,12 +110,12 @@ public void testW_trig() {
  .startPosition(1)
  .trigonometric(53)
  .build();
- for (int i = 0; i < 25; ++i) {
+// for (int i = 0; i < 25; ++i) {
 // System.out.println(ch.test(1+2 * i, 2));
 // }
 // System.out.println(ch.testAll());
 // System.out.println(ch.robustTestCoefficients(i*2, 2).getValue());
- }
+// }
  }
 
  @Test

jdplus-main-base/jdplus-sa-base-parent/jdplus-sa-base-r/src/main/java/jdplus/sa/base/r/SeasonalityTests.java

@@ -122,26 +122,33 @@ public double[] canovaHansen(double[] s, int period, boolean trig, boolean lag1,
  .startPosition(start);
  if (trig) {
  CanovaHansen ch = builder.trigonometric(period).build();
- double[] q = new double[1 + period / 2];
  boolean even = period % 2 == 0;
- int nq = even ? q.length - 2 : q.length - 1;
+ int p2=period/2;
+ int nq = even ? p2-1 : p2;
+ double[] q = new double[3 + p2];
  int icur = 0;
  for (int i = 0; i < nq; ++i, ++icur) {
  q[icur] = ch.test(i * 2, 2);
  }
  if (even) {
  q[icur++] = ch.test(period - 2);
  }
- q[icur] = ch.testAll();
+ q[icur++] = ch.testAll();
+ StatisticalTest seasonalityTest = ch.seasonalityTest();
+ q[icur++]=seasonalityTest.getValue();
+ q[icur]=seasonalityTest.getPvalue();
  return q;
  } else {
  CanovaHansen ch = builder.dummies(period).build();
- double[] q = new double[period + 1];
+ double[] q = new double[period + 3];
  for (int i = 0; i < period; ++i) {
  q[i] = ch.test(i);
  }
  q[period] = ch.testAll();
- return q;
+ StatisticalTest seasonalityTest = ch.seasonalityTest();
+ q[period+1]=seasonalityTest.getValue();
+ q[period+2]=seasonalityTest.getPvalue();
+ return q;
  }
  }
 

jdplus-main-base/jdplus-sa-base-parent/jdplus-sa-base-r/src/test/java/jdplus/sa/base/r/SeasonalityTestsTest.java

@@ -72,10 +72,10 @@ public void testCombinedTest() {
 
  @Test
  public void testCH() {
- DoubleSeq x=DoubleSeq.of(Data.ABS_RETAIL).log();
- double[] q=SeasonalityTests.canovaHansen(x.toArray(), 12, false, true, "Bartlett", 15, 3);
-// System.out.println(DoubleSeq.of(q));
- q=SeasonalityTests.canovaHansen(x.toArray(), 12, true, true, "Bartlett", 15, 3);
-// System.out.println(DoubleSeq.of(q));
+ DoubleSeq x=DoubleSeq.of(Data.RETAIL_BOOKSTORES).log().delta(1);
+ double[] q=SeasonalityTests.canovaHansen(x.toArray(), 12, false, false, "Bartlett", -1, 3);
+ System.out.println(DoubleSeq.of(q));
+ q=SeasonalityTests.canovaHansen(x.toArray(), 12, true, false, "Bartlett", -1, 3);
+ System.out.println(DoubleSeq.of(q));
  }
 }

jdplus-main-base/jdplus-toolkit-base-parent/jdplus-toolkit-base-core/src/main/java/jdplus/toolkit/base/core/data/analysis/WindowFunction.java

@@ -36,23 +36,27 @@ public enum WindowFunction {
 
  /**
  * Returns the normalized window function, defined on [-1, 1].
- * We must have that f(-1)=f(1)=0
+ * We must have that f(-1)=f(1)=0, f(0)=1
  * The window function is even, so that f(-x)=f(x)
  * It should be noted that the functions don't check the validity of the input
  *
  * @return
  */
  public DoubleUnaryOperator window() {
  switch (this) {
- case Welch:
+ case Welch -> {
  return x -> 1.0 - x * x;
- case Tukey:
+ }
+ case Tukey -> {
  return x -> 0.5 * (1 + Math.cos(Math.PI * x));
- case Bartlett:
+ }
+ case Bartlett -> {
  return x -> x < 0 ? 1 + x : 1 - x;
- case Hamming:
+ }
+ case Hamming -> {
  return x -> 0.54 + 0.46 * Math.cos(Math.PI * x);
- case Parzen:
+ }
+ case Parzen -> {
  return x -> {
  double x1 = x < 0 ? -x : x;
  if (x <= .5) {
@@ -63,8 +67,10 @@ public DoubleUnaryOperator window() {
  return 2.0 * y * y * y;
  }
  };
- case Square:
+ }
+ case Square -> {
  return x -> 1;
+ }
  }
  return null;
  }