/* MSC-P-033 marketing forecast validation. MIT License. Secondary implementation.
   The design matrix, the origins, the horizon and the nominal level are the
   declared ones, but the least-squares fit is delegated to PROC REG, whose
   solver and residual conventions are its own, so the last decimals may differ
   from the Python, R and SPSS references. This file is an independent check of
   the two flags and of the verdict, not a digit-for-digit reproduction. */
%let period = 52;
%let harmonics = 2;
%let first_origin = 104;
%let origin_step = 4;
%let last_origin = 152;
%let horizon = 4;
%let nominal = 0.80;
%let normal_quantile = 1.281552;  /* the 0.90 quantile of the standard normal law */
%let mase_max = 1.00;
%let coverage_tolerance = 0.10;

filename p033csv "public/datasets/msc-p033-weekly-series.csv";
data _null_;
  infile p033csv obs=1 lrecl=32767 truncover;
  input;
  if strip(_infile_) ne "week,revenue_eur" then do;
    put "ERROR: exact schema required"; abort cancel;
  end;
run;
data p033;
  infile p033csv dsd firstobs=2 truncover lrecl=32767 end=eof;
  input @;
  if countc(_infile_, ',') ne 1 then do; put "ERROR: exactly two cells required per row"; abort cancel; end;
  input week revenue_eur;
  if week ne _n_ then do; put "ERROR: weeks must be numbered 1, 2, ... without gaps"; abort cancel; end;
  if missing(revenue_eur) or revenue_eur <= 0 then do; put "ERROR: the outcome must be finite and strictly positive"; abort cancel; end;
  /* Declared design: intercept, linear trend and two annual harmonics. */
  trend = week;
  sin1 = sin(2 * constant('pi') * 1 * week / &period);
  cos1 = cos(2 * constant('pi') * 1 * week / &period);
  sin2 = sin(2 * constant('pi') * 2 * week / &period);
  cos2 = cos(2 * constant('pi') * 2 * week / &period);
  if eof then call symputx("weeks", _n_);
run;
%if &weeks < %eval(&last_origin + &horizon) %then %do;
  %put ERROR: the series is shorter than the declared validation design; %abort cancel;
%end;

%macro validate;
  data evaluations; length flag $4; stop; run;
  %local origin h;
  %do origin = &first_origin %to &last_origin %by &origin_step;
    /* Nothing beyond the origin is read by the fit. */
    data window; set p033; if week <= &origin; run;
    proc reg data=window noprint outest=estimates edf;
      model revenue_eur = trend sin1 cos1 sin2 cos2;
    run; quit;
    data _null_;
      set estimates;
      call symputx("b0", intercept);
      call symputx("b1", trend);
      call symputx("b2", sin1);
      call symputx("b3", cos1);
      call symputx("b4", sin2);
      call symputx("b5", cos2);
      call symputx("sigma", _rmse_);
    run;
    /* Declared scale: in-sample mean absolute seasonal naive error. */
    proc sql noprint;
      select mean(abs(a.revenue_eur - b.revenue_eur)) into :scale
      from p033 a inner join p033 b on a.week = b.week + &period
      where a.week <= &origin;
    quit;
    data step;
      length flag $4;
      set p033;
      if &origin < week <= &origin + &horizon;
      horizon = week - &origin;
      prediction = &b0 + &b1 * trend + &b2 * sin1 + &b3 * cos1 + &b4 * sin2 + &b5 * cos2;
      absolute = abs(revenue_eur - prediction);
      scaled = absolute / &scale;
      half_width = &normal_quantile * &sigma;
      lower = prediction - half_width;
      upper = prediction + half_width;
      inside = (lower <= revenue_eur <= upper);
      width = upper - lower;
      score = width
            + (2 / (1 - &nominal)) * max(0, lower - revenue_eur)
            + (2 / (1 - &nominal)) * max(0, revenue_eur - upper);
      keep week horizon absolute scaled inside width score;
    run;
    /* The benchmark is the seasonal naive forecast on the same scale. */
    proc sql noprint;
      create table benchmark_step as
        select a.week, abs(a.revenue_eur - b.revenue_eur) / &scale as scaled_benchmark
        from p033 a inner join p033 b on a.week = b.week + &period
        where &origin < a.week <= &origin + &horizon;
    quit;
    data step; merge step benchmark_step; by week; run;
    proc append base=evaluations data=step force; run;
  %end;
  proc sql noprint;
    select mean(scaled), mean(scaled_benchmark), mean(absolute), mean(inside), mean(width), mean(score)
      into :mase, :benchmark, :mae, :coverage, :width, :score
    from evaluations;
  quit;
  proc means data=evaluations noprint nway;
    class horizon; var scaled;
    output out=by_horizon(drop=_type_ _freq_) mean=mase_by_horizon;
  run;
  data summary;
    length accuracy_flag calibration_flag $4 verdict $34;
    forecasts = &sqlobs;
    model_mase = &mase;
    benchmark_mase = &benchmark;
    model_mae = &mae;
    empirical_coverage = &coverage;
    mean_interval_width = &width;
    mean_interval_score = &score;
    accuracy_flag = ifc(model_mase < &mase_max, "PASS", "FAIL");
    calibration_flag = ifc(abs(empirical_coverage - &nominal) <= &coverage_tolerance, "PASS", "FAIL");
    verdict = ifc(accuracy_flag = "PASS" and calibration_flag = "PASS",
                  "FORECAST_READABLE_FOR_PLANNING", "DIAGNOSTIC_BLOCKS_FORECAST_READING");
  run;
  proc print data=summary noobs; run;
  proc print data=by_horizon noobs; run;
%mend validate;
%validate;