*
* @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/2005 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[real] w
local integer i q d startl endl
local real avgw nonpos
*
inquire(seasonal) d<>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