*
* @ClassicalDecomp(options)  x  start  end
*  Does a classical decomposition of the series "x" into
*  a trend-cycle, seasonal and irregular. The seasonals are
*  assumed to be constant across the data range.
*
*  With SEASONAL=ADDITIVE (the default), the seasonal factors sum
*  to zero, and the data are assumed to take the form
*
*   X=TC+S+I
*
*  With SEASONAL=MULTIPLICATIVE, the seasonal factors are scaled to
*  average 1.0 and the data are assumed to take the form
*
*   X=TCxSxI
*
*  The TREND option determines how the final trend-cycle is estimated
*  after the seasonal factors have been computed. LINEAR and QUADRATIC
*  request regressions with 1, t or 1, t and t**2 respectively. The
*  default (TREND=MA) does a moving average with correction for end
*  terms.
*
* Options:
*   SEASONAL=[ADDITIVE]/MULTIPLICATIVE
*   TREND=NONE/LINEAR/QUADRATIC/[MA]
*   SPAN=seasonal span [CALENDAR seasonal]
*   TC=estimated trend cycle
*   FACTORS=estimated seasonal factors
*   IRREG=estimated irregular component
*
* Revision Schedule:
*  06/05 Written by Tom Doan, Estima
*
proc ClassicalDecomp x start end
type series x
type integer start end
*
option integer span
option choice  trend     3  none linear quadratic ma
option choice  seasonal  1  additive multiplicative
option series  *factors
option series  *tc
option series  *irreg
*
local series rawseason deseason tcinit trend1 trend2
local vect[integer] lags
local vect[real] coeffs
local vect[real] w
local integer i q d startl endl
local real avgw nonpos
*
inquire(seasonal) d<<span
if d==0 {
   disp 'Syntax: @ClassicalDecomp(span=seasonal) input start end output'
   return
}

inquire(series=x) startl<<start endl<<end
sstats startl endl (x<=0)>>nonpos
if seasonal==2.and.nonpos>0
  disp "@ClassicalDecomp(seasonal=mult) can't be used with non-positive data"
*
* Compute a centered, flat MA of width "d" periods to estimate
* the trend cycle
*
compute q=d/2
filter(type=centered,width=d) x startl+q endl-q tcinit
*
* Determine raw seasonal
*
if seasonal==1
   set rawseason startl+q endl-q = x-tcinit
else
   set rawseason startl+q endl-q = x/tcinit
*
* Compute deviations from TC for each subperiod
*
dim w(d)
do i=1,d
   sstats(mean,smpl=%clock(t,d)==i) startl+q endl-q rawseason>>w(i)
end do i
*
* Center the seasonals
*
if seasonal==1 {
   compute avgw=%avg(w)
   compute w=w-avgw
   set deseason startl endl = x-w(%clock(t,d))
}
else {
   compute avgw=%avg(w)
   compute w=w/avgw
   set deseason startl endl = x/w(%clock(t,d))
}
*
* Extract the seasonal
*
if %defined(factors)
   set factors startl endl = w(%clock(t,d))

if trend==4 {
   *
   * Do an MA with asymmetrical window to cover the end terms
   *
   graph
   # x
   filter(type=centered,width=d,extend=mse) deseason startl endl tc
   if %defined(irreg)
      set irreg startl endl = x-tc
}
else {
   *
   * Regress on requested trend variables, returning the residual.
   *
   set trend1 startl endl   = t
   set trend2 startl endl   = t**2
   linreg(entries=trend,noprint) deseason startl endl irreg
   # constant trend1 trend2
   if %defined(tc)
      prj tc startl endl
}
end


