*
* @CFeat( options )  series  start  end
*
* Program for the determination of the turning points of the time series,
* calculation of the duration and amplitude of the expansions and recessions.
* The expansion (recession) is defined as a timespan between a cyclical
* trough (peak) and peak (trough). The amplitude hence is the absolute distance
* between trough (peak) and peak (trough).
*
* The output of the program is the table of the durations and amplitudes of the
* cyclical phases of the particular cycles as well as the average durations and
* amplitudes.
*
* Parameters:
*   series      input series
*   start end   range to analyze
*
* Options:
*    m = minimum cycle duration [8]
*    r = determines (inversely) the minimum magnitude of an extremum [1]
*    asym = output vector (dimension 6), containing the average durations
*          and amplitudes (in pairs of duration and amplitude) of the
*          whole cycle, expansion, and recession.
*    table=none/[full]/summary
*
* Revision Schedule:
*   05/2000 Written by K.A.Kholodilin, UAB
*   06/2002 Revised by K.A.Kholodilin, UAB
*   04/2005 Revised by Tom Doan, Estima.
*      Syntax changed to use options rather than parameters, REPORT instruction
*      replaces DISPLAY, consecutive peaks or troughs with no intervening opposite
*      cycle merged.
*
procedure cfeat series startp endp
type series series
type integer startp endp
*
option integer m    8
option integer r    1
option vector  *asym
option choice  table 2 none full summary
*
local vec x
local vec chron chro
local rec tabl dur
local series points
local integer i j d k n np nexp nrec
local real durexp durrec amexp amrec avdexp avdrec avamexp avamrec
local integer start end start1 end1 count cnt
local real stdv
local integer lastkept
local real lastval

inquire(series=series) start<<startp end<<endp
comp n=end-start
comp k=n-m+1
comp d=fix((m-1)/2)

dim x(m) chron(n) chro(n)

comp chron = %const(0)
statistics(noprint) series
comp stdv=%variance**0.5

do i=1,k
   do j=1,m
      comp x(j)=series(i+j-1)
   end do j
   if x(d)>=%maxvalue(x)
      comp chron(i)=x(d)
   else if x(d)<=%minvalue(x)
      comp chron(i)=x(d)
   else
      comp chron(i)=0.0
end do i
*
* Check for peaks and troughs whose deviation is smaller than
* (1/r) x the standard deviation of the original series. Also, merge
* consecutive peaks or troughs (for instance, two local peaks without an
* intervening cyclical trough), keeping the more extreme of the two.
*
compute lastkept=0,lastval=0.0
do i=1,k
   if abs(chron(i))>=stdv/r {
      comp chron(i)=chron(i)
      if lastval*chron(i)>0 {
         if abs(chron(i))<abs(lastval) {
            compute chron(i)=0.0
            next
         }
         else
            compute chron(lastkept)=0.0
      }
      compute lastval=chron(i),lastkept=i
   }
   else
      comp chron(i)=0.0
end do i
*
* Redate the information
*
set points = 0.0
comp start1=start+d-1
do i=1,k
   comp end1=start1
   set points start1 end1 = chron(i)
   comp start1=start1+1
end do i
*
* Count the number of cycles to dimension the table matrix
*
comp count=0
do i=1,k
   if abs(chron(i))>0
      comp count=count+1
   else
      comp count=count
end do i

dim tabl(count,4)
comp tabl=%const(0)
comp cnt=0

do i=1,k
   if chron(i)>0
      comp cnt=cnt+1, tabl(cnt,1)=i+d-1, tabl(cnt,2)=chron(i)
   else if chron(i)<0
      comp cnt=cnt+1, tabl(cnt,3)=i+d-1, tabl(cnt,4)=chron(i)
   else
      comp cnt=cnt
end do i

do j=1,cnt-1
   if tabl(j,2)>tabl(j+1,2)
      comp tabl(j,2)=tabl(j,2), tabl(j+1,2)=0
   else
      comp tabl(j,2)=0, tabl(j,1)=0
   if tabl(j,4)<tabl(j+1,4)
      comp tabl(j,4)=tabl(j,4), tabl(j+1,4)=0
   else
      comp tabl(j,4)=0, tabl(j,3)=0
end do j

comp np=fix(cnt/2)
dim dur(np,8)
comp dur=%const(0)

do i=1,cnt-1
   comp k=fix((i+1)/2)
   if tabl(i,1)>0
      comp dur(k,5)=tabl(i,1), dur(k,6)=tabl(i+1,3), dur(k,7)=tabl(i+1,3)-tabl(i,1), dur(k,8)=tabl(i,2)-tabl(i+1,4)
   else
      comp dur(k,1)=tabl(i,3), dur(k,2)=tabl(i+1,1), dur(k,3)=tabl(i+1,1)-tabl(i,3), dur(k,4)=tabl(i+1,2)-tabl(i,4)
end do i

comp nexp=0, nrec=0, durexp=0.0, durrec=0.0, amexp=0.0, amrec=0.0
do i=1,np
   if dur(i,1)>0
      comp nexp=nexp+1, durexp=durexp+dur(i,3), amexp=amexp+dur(i,4)
   if dur(i,5)>0
      comp nrec=nrec+1, durrec=durrec+dur(i,7), amrec=amrec+dur(i,8)
end do i

comp avdexp=durexp/nexp, avdrec=durrec/nrec, avamexp=amexp/nexp, avamrec=amrec/nrec

if table==2
{
   report(action=define)
   report(row=new,atcol=1,span) "Business Cycle Phases of "+%l(series)
   report(row=new,atcol=1) "Expansions" "" "" "" "Recessions"
   report(row=new,atcol=1) "Start" "End" "Duration" "Amplitude" "Start" "End" "Duration" "Amplitude"
   do i=1,np
      report(row=new)
      if dur(i,1)>0
         report(row=current,atcol=1) %datelabel(fix(dur(i,1))) %datelabel(fix(dur(i,2))) dur(i,3) dur(i,4)
      if dur(i,5)>0
         report(row=current,atcol=5) %datelabel(fix(dur(i,5))) %datelabel(fix(dur(i,6))) dur(i,7) dur(i,8)
   end do i
   report(row=new,atcol=1) "Average duration"
   report(row=current,atcol=3) avdexp "" "" "" avdrec
   report(row=new,atcol=1) "Average amplitude"
   report(row=current,atcol=4) avamexp "" "" "" avamrec
   report(action=format,atrow=3,atcol=3,tocol=3,picture="*.##")
   report(action=format,atrow=3,atcol=7,tocol=7,picture="*.##")
   report(action=format,atrow=3,atcol=4,tocol=4,width=8)
   report(action=format,atrow=3,atcol=8,tocol=8,width=8)
   report(action=show)
}

else
if table==3
{
   display  @30 "Business Cycle Phases" %label(series)
   display @10 "Expansions" @50 "Recessions"
   display @10 "Average duration" @30 avdexp @50 avdrec
   display @10 "Average amplitude" @40 avamexp @60 avamrec
}

if %defined(asym) {
   dim asym(6)
   comp asym(1)=avdexp+avdrec
   comp asym(2)=avamexp+avamrec
   comp asym(3)=avdexp
   comp asym(4)=avamexp
   comp asym(5)=avdrec
   comp asym(6)=avamrec
}

end cfeat




