ocsb function

OCSBResult ocsb(List<double> x, int s, { int? kLags, })

Implementation

OCSBResult ocsb(List<double> x, int s, {int? kLags}){
  final n=x.length; if (n < s+10) return OCSBResult(double.nan, 0);
  // following Hyndman-Khandakar simplified OCSB:
  // ΔX_t^s = β X_{t-s} + Σ φ_i ΔX_{t-i} + ε_t
  final dx = differenceSeas(x, s, 1);
  final d1 = differenceOrd(x, 1);
  final m = dx.length;
  final k = kLags ?? math.max(0, (math.pow(n/100.0, 0.25)).floor()).toInt();
  final rows = m - math.max(k,1).toInt();
  if (rows <= 3) return OCSBResult(double.nan, k);
  final X = List<List<double>>.generate(rows, (_)=> <double>[]);
  final Y = List<double>.filled(rows, 0.0);
  for (int t=math.max(k,1); t<m; t++){
    final row=<double>[1.0, dx[t-1]]; // const + X_{t-1}^s
    for (int i=1;i<=k;i++){
      final di = (t-i>=0 && t-i < d1.length) ? d1[t-i] : 0.0;
      row.add(di);
    }
    X[t-math.max(k,1)] = row;
    Y[t-math.max(k,1)] = dx[t];
  }
  final beta = _solveDense( // beta = [c, β, φ1..]
    (){ final kcols=X[0].length; final XtX=List.generate(kcols,(_)=>List<double>.filled(kcols,0.0));
      final XtY=List<double>.filled(kcols,0.0);
      for (int i=0;i<rows;i++){ for (int a=0;a<kcols;a++){ XtY[a]+=X[i][a]*Y[i]; for (int b=0;b<kcols;b++) XtX[a][b]+=X[i][a]*X[i][b]; } }
      return XtX; }(),
    (){ final kcols=X[0].length; final XtY=List<double>.filled(kcols,0.0);
      for (int i=0;i<rows;i++) for (int a=0;a<kcols;a++) XtY[a]+=X[i][a]*Y[i];
      return XtY; }()
  );
  // get standard error of β (index 1)
  final kcols=X[0].length;
  final yhat=List<double>.generate(rows,(i){ var s=0.0; for (int j=0;j<kcols;j++) s+=X[i][j]*beta[j]; return s; });
  var sse=0.0; for (int i=0;i<rows;i++) sse+=(Y[i]-yhat[i])*(Y[i]-yhat[i]);
  final sigma2 = sse / (rows - kcols);
  // invert XtX
  final inv = (){
    final XtX=List.generate(kcols,(_)=>List<double>.filled(kcols,0.0));
    for (int i=0;i<rows;i++) for (int a=0;a<kcols;a++) for (int b=0;b<kcols;b++) XtX[a][b]+=X[i][a]*X[i][b];
    // Gauss-Jordan
    final A=List.generate(kcols,(i)=>List<double>.from(XtX[i])); final I=List.generate(kcols,(i)=>List<double>.filled(kcols,0.0)); for (int i=0;i<kcols;i++) I[i][i]=1.0;
    for (int i=0;i<kcols;i++){
      int piv=i; double mx=A[i][i].abs(); for (int r=i+1;r<kcols;r++){ final v=A[r][i].abs(); if (v>mx){mx=v; piv=r;} }
      final t=A[i]; A[i]=A[piv]; A[piv]=t; final ti=I[i]; I[i]=I[piv]; I[piv]=ti;
      final d=A[i][i]; for (int c=0;c<kcols;c++){ A[i][c]/=d; I[i][c]/=d; }
      for (int r=0;r<kcols;r++){ if (r==i) continue; final f=A[r][i]; for (int c=0;c<kcols;c++){ A[r][c]-=f*A[i][c]; I[r][c]-=f*I[i][c]; } }
    }
    return I;
  }();
  final seBeta = math.sqrt(sigma2 * inv[1][1]);
  final tstat = beta[1] / (seBeta==0?1e-12:seBeta);
  return OCSBResult(tstat, k);
}