/***********************************************************************************************/
/*                                                                                             */
/*  Modified BBQ program                                                                       */ 
/*                                                                                             */ 
/*  Computes turning points and imposes restrictions in one step                               */ 
/*                                                                                             */ 
/*    NB: rule is peak greater than turnphase on either side, trough less than turnphase       */
/*				 on either side.             (can be modified)                                       */ 
/*                                                                                             */ 
/*    Restriction: Phase quarters/months, Cycle  quarters/months Alternating Troughs and Peaks,*/ 
/*                  if two peaks(troughs) in a row chooses highest(lowest), Trough must be     */ 
/*                  lower then preceeding peak, and No truning point within phase length       */ 
/*                  of end points    + others that are needed    						              	   */
/*																						                                            	   */
/*                                                                                             */
/*    Date: 22th March 2005                                                                    */
/*                                                                                             */
/*    Author: James Engel - Code modified from Adrian Pagan and Don Harding's BBQ code         */
/*                                                                                             */ 
/***********************************************************************************************/
cls;
library pgraph;
ts = hsec; 

nrep= 1;       @ set to one if analysing real data  - line 105 to enter data/model @
complete = 0;  @ switch: 1-only use complete cycles,0-use incomplete cycles (excess still on complete cycle)@
                   
                               
@ Cycle Characteristics @

turnphase = 2;      @ initial turning point peak/trough decision @
phase = 2;          @ censoring rules @
cycle = 5;
thresh = 10.4;       @ bypasses phase and cycle restriction if peak to trough is > than thresh @

@ Frequency @
freq = 1;       @ 1 for quarterly, 2 for monthly @
                @ note - enter data or model line 110 or so @


ndfy = 1947;             @ First Year of data set @
ndfm = 1;                @ First quarter/month of data set @
ndly = 2002;             @ Last Year of data set @
ndlm = 2;                @ Last quarter/month of data set @


/*ndfy = 1970;             @ First Year of data set @
ndfm = 1;                @ First quarter/month of data set @
ndly = 2003;             @ Last Year of data set @
ndlm = 4;                @ Last quarter/month of data set @*/




nsfy=ndfy;nsfm=ndfm;nsly=ndly;nslm=ndlm;

if freq ==1;
nd= 4*(ndly-1-ndfy) + (5-ndfm) + (ndlm);     @ Number of data points @
elseif freq == 2;
nd= 12*(ndly-1-ndfy) + (13-ndfm) + (ndlm); 
endif;

nn=nd;


notentp=0;
durs=zeros(nd,1);bcp5=zeros(nd,1);bct5=zeros(nd,1);st=zeros(nd,1);

pds=zeros(1,1);tds=zeros(1,1);pdsa=zeros(1,1);tdsa=zeros(1,1);bvec=zeros(nd,1);

pdmmat=zeros(nrep+1,1);tdmmat=zeros(nrep+1,1);pdmmata=zeros(nrep,1);tdmmata=zeros(nrep,1);

@set icensor@
@icensor =0 if no censoring =1 if is@
icensor=1;
@initialize random numbers@

rndseed 7654;



PDCV  =0; TDCV  =0;PACV =0;TACV =0;PDM =0;PDMA=0; TDM =0; TDMA  =0;PDCM=0; TDCM=0;
PDEM =0;    TDEM =0;    TDEMA=0;    PDEMA =0; EPCV = 0; ETCV = 0;

nbt=0;
nbp=0;
iter=1;


"NO ITERS";NREP;
do while iter<=nrep;

nd=nn;

/* Data Input - needs to be in LN form */

@load x[nd,1]=ozgdp.txt;@
@load x[nd,1]=usgdp.txt;@
@load x[nd,1]=cangdp.txt;@
@load x[nd,1]=c:\rotterdam\eurogdp.dat;@

@x = ln(x);@


/* Simulated Model data - place model at end of code */

@{x}=usar;@




y=zeros(nd,1);c=zeros(nd,1);in=zeros(nd,1);
nbp=0;nbt=0;
{bcp5,bct5,nbp,nbt}=rawall(x[1:nd]);  @ calculates turning points with restrictions @

if nbp+nbt<=2;
notentp=notentp+1;
goto skipup;
endif;
update(x,bcp5,bct5,nd,nbp,nbt);
skipup:

if  nrep == 1.and(nbp+nbt>2);
@nbp;nbt;@  @ number of peaks and number of troughs @

"peaks at";
bcp5[1:nbp];

"troughs at";
bct5[1:nbt];
?;

@histo(x,bcp5,bct5,nd,nbp,nbt); @   @ to graph expansions @
?;
endif;

pdmmat[iter+1]=pdm-pdmmat[iter];
tdmmat[iter+1]=tdm-tdmmat[iter];
iter=iter+1;

endo;



if nrep == 1.and(nbp+nbt>2);
/*remove below if want states */
st=zeros(nd,1);


{st}=states(bcp5,bct5,nbp,nbt);

@determine obs in which states have been completed@
na=minc(bct5[1]|bcp5[1]);
nb=maxc(bct5[nbt]|bcp5[nbp]);
nb-na+1;
z=x[na:nb]~st[na:nb];
z1 = seqa(ndfy+(ndfm-1)/4,0.25,nd)~exp(x)~st;
@"";
"Dated Series"
"";@
z1;

@output file=c:\unsw2005\econ5176\cstates.dat reset;@
output file=c:\rotterdam\infstates.dat reset;
output on;
st;
output off;
@ Writes data and states to spreadsheet @

@ret = xlswritem(z1,"C:\unsw2005\econ5176\StatesData" ,"a1", 1, 0); @

endif;



nrep1=nrep-notentp;
?;
"statistics on average cycle";
"contractions/expansions";

"durations";

pdm/nrep1;tdm/nrep1;

"amplitudes";

pdma/nrep1;tdma/nrep1;

"cumulative movements";

pdcm/nrep1;tdcm/nrep1;

"excess movements % of triange area";

pdem/nrep1;tdem/nrep1;

"cv of dur";
pdcv/nrep1;tdcv/nrep1;
"cv of amps";
pacv/nrep1;tacv/nrep1;
"cv of xss";
epcv/nrep1;etcv/nrep1;


"no of its skipped since no peaks+troughs <=2";notentp;
pa=sortc(pdmmat,1);
pb=sortc(tdmmat,1);

et = hsec-ts;
"running time";
et/100;

screen off;
pb;


end;




/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/




proc(0) = update(x,bcp5,bct5,nd,nbp,nbt);


local ntr, npk;
local td,pd,pdc,tdc,pde,tde,pda,tda, tdea, pdea;
local za,r,k,nr,nv,IR,r1;






@ p stands for peaks, t for troughs, p gives cntractions, t expansions@

@code is pd,td is durations; pda,tda is amps; pdc,tdc is cum move pde,
 tde exces@

@excess is measured differently to avoid case that amps is close to zero
in part cycle so that denom can become neg. So use cum movements as denom
now vs triangle in early paper@

@ Calculate Peak To Trough Durations and amps @
@ p stands for peaks, t for troughs, p gives cntractions, t expansions@

@code is pd,td is durations; pda,tda is amps; pdc,tdc is cum move pde,
 tde exces@

ntr=nbt;npk=nbp;

nr=nbt|nbp;
nv=maxc(nr);

pdc=zeros(nv,1);tdc=zeros(nv,1);
PDA=ZEROS(Nv,1);tda=zeros(nv,1);td=zeros(nv,1);pd=zeros(nv,1);
pdc=zeros(nv,1);tdc=zeros(nv,1);pde=zeros(nv,1);tde=zeros(nv,1);
tdea=zeros(nv,1);pdea=zeros(nv,1);


@ Calculate Peak To Trough Durations and amps @
if bcp5[1,1] .< bct5[1,1]; @ Peaks are first @

    nr=nbt|nbp;
    r=NBT;
     pd=bct5[1:r,1]-bcp5[1:r,1];
     pda=x[bct5[1:r,1]]-x[bcp5[1:r,1]];

    k=1;
    do while k<=r;
    
    pdc[k]=sumc(x[bcp5[k,1]:bct5[k,1],1]-x[bcp5[k,1],1]);
    k=k+1;
    endo;

else;                      @ troughs are First@

    r=NBT-1;
    pd=bct5[2:r+1,1]-bcp5[1:r,1];
    pda=x[bct5[2:r+1,1]]-x[bcp5[1:r,1]];
    
    k=1;
    do while k<=r;
     pdc[k]=sumc(x[bcp5[k,1]:bct5[k+1,1],1]-x[bcp5[k,1],1]);
    
    k=k+1;
    endo;
    
    r1=r;

endif;


@ Calculate Trough to Peak Durations and amplitudes@
if bct5[1,1] .< bcp5[1,1]; @ Troughs are first @

    r=NBP;
     td=bcp5[1:r,1]-bct5[1:r,1];
     tda=x[bcp5[1:r,1]]-x[bct5[1:r,1]];
    
    k=1;
    do while k<=r;
     tdc[k]=sumc(x[bct5[k,1]:bcp5[k,1],1]-x[bct5[k,1],1]);
    
    k=k+1;
    endo;

else;                      @ peaks are First @

     r=NBP-1;
    
     td=bcp5[2:r+1,1]-bct5[1:r,1];
     tda=x[bcp5[2:r+1,1]]-x[bct5[1:r,1]];
    
    k=1;
    do while k<=r;
        
     tdc[k]=sumc(x[bct5[k,1]:bcp5[k+1,1],1]-x[bct5[k,1],1]);
        
    k=k+1;
    endo;


endif;



pdc=pdc[1:rows(pd)];
tdc=tdc[1:rows(td)];

@compute excesses note it is now a percentage of triangle area@

za=(pd.*pda)/2;
pde=(pdc-za-.5*pda)./za;
pdea=100*abs(pdc-((pd.*pda)/2))./pd;
pdea=100*(pdc-((pd.*pda)/2))./za;

za=(td.*tda)/2;
tde=(tdc-za-.5*tda)./za;
tdea=100*abs(tdc-((td.*tda)/2))./td;
tdea=100*(tdc-((td.*tda)/2))./za;


/*****************************************************************/

if complete == 0;    @ switch: 1-only use complete cycles,0-use incomplete cycles (excess still on complete cycle)@

        if bcp5[1,1] .< bct5[1,1];     @ modifies code to include incomplete cycles @
        
        	bct5  = 1|bct5;
            goto endb;
        elseif 	bct5[1,1] .< bcp5[1,1];
        
        	bcp5 = 1|bcp5; 
            goto endb;
        endif;
        
        endb:
        nbt = rows(bct5);
        nbp = rows(bcp5);
        
        if bcp5[nbp,1] .< bct5[nbt,1];     @ modifies code to include incomplete cycles @
        
        	bcp5  = bcp5|nd;
            goto ende;
        elseif 	bct5[nbt,1] .< bcp5[nbp,1];
        
        	bct5 = bct5|nd; 
            goto ende;
        endif;
        
        ende:
        
        nbt = rows(bct5);
        nbp = rows(bcp5);

        ntr=nbt;npk=nbp;
        
        nr=ntr|npk;
        nv=maxc(nr);
        
        pdc=zeros(nv,1);tdc=zeros(nv,1);
        PDA=ZEROS(Nv,1);tda=zeros(nv,1);td=zeros(nv,1);pd=zeros(nv,1);
        pdc=zeros(nv,1);tdc=zeros(nv,1);
        
        @ Calculate Peak To Trough Durations and amps @
        if bcp5[1,1] .< bct5[1,1]; @ Peaks are first @
        
            nr=ntr|npk;
            r=ntr;
             pd=bct5[1:r,1]-bcp5[1:r,1];
             pda=x[bct5[1:r,1]]-x[bcp5[1:r,1]];
        
            k=1;
            do while k<=r;
            
            pdc[k]=sumc(x[bcp5[k,1]:bct5[k,1],1]-x[bcp5[k,1],1]);
            k=k+1;
            endo;
        
        else;                      @ troughs are First@
        
            r=ntr-1;
            pd=bct5[2:r+1,1]-bcp5[1:r,1];
            pda=x[bct5[2:r+1,1]]-x[bcp5[1:r,1]];
            
            k=1;
            do while k<=r;
             pdc[k]=sumc(x[bcp5[k,1]:bct5[k+1,1],1]-x[bcp5[k,1],1]);
            
            k=k+1;
            endo;
            
            r1=r;
        
        endif;
        
        
        @ Calculate Trough to Peak Durations and amplitudes@
        if bct5[1,1] .< bcp5[1,1]; @ Troughs are first @
        
            r=npk;
             td=bcp5[1:r,1]-bct5[1:r,1];
             tda=x[bcp5[1:r,1]]-x[bct5[1:r,1]];
            
            k=1;
            do while k<=r;
             tdc[k]=sumc(x[bct5[k,1]:bcp5[k,1],1]-x[bct5[k,1],1]);
            
            k=k+1;
            endo;
        
        else;                      @ peaks are First @
        
             r=npk-1;
            
             td=bcp5[2:r+1,1]-bct5[1:r,1];
             tda=x[bcp5[2:r+1,1]]-x[bct5[1:r,1]];
            
            k=1;
            do while k<=r;
                
             tdc[k]=sumc(x[bct5[k,1]:bcp5[k+1,1],1]-x[bct5[k,1],1]);
                
            k=k+1;
            endo;
        
        
        endif;
        
        
        
        pdc=pdc[1:rows(pd)];
        tdc=tdc[1:rows(td)];

endif;
/**************************************************************/

@cumulate. When nrep=1 then gives raw data, otherwsise sums over monte
carlo amounts@

 
pdcv=pdcv+stdc(pd)/meanc(pd);       @compute cv's@
tdcv=tdcv+stdc(td)/meanc(td);
pacv=pacv+stdc(pda)/meanc(pda);
tacv=tacv+stdc(tda)/meanc(tda);

epcv = epcv + stdc(pde)/meanc(pde); @ cv of xss @
etcv = etcv + stdc(tde)/meanc(tde);




pdm=pdm+meanc(pd);          @ durations @
pdma=pdma+meanc(pda);

tdm=tdm+meanc(td);          @ durations @
tdma=tdma+meanc(tda);


pdcm=pdcm+meanc(pdc);       @ cumulative  @
tdcm=tdcm+meanc(tdc);


tdem=tdem+meanc(tde);       @ xss @
tdema=tdema+meanc(tdea);
pdem=pdem+meanc(pde);
pdema=pdema+meanc(pdea);

retp;
endp;


/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/


proc(1) = states(bcp5,bct5,nbp,nbt);
local nrp,i,nrt, nrpp;

@1=expansion;0=contraction@


if bcp5[1] < bct5[1];
goto touch;
endif;
i=1;
nrp=nbp;
do while i<=nrp;
@  bct5[i];bcp5[i];@
st[bct5[i]+1:bcp5[i]]=ones(bcp5[i]-bct5[i],1);
i=i+1;
endo;

@test if one more troughs than peaks if so then pattern is
tptpt and need to fill in with ones at end@



if nbt>nbp;
st[bct5[nbp+1]+1:nd]=ones(nd-bct5[nbp+1],1);
endif;

goto ret;
touch:
i=2;

@need to determine if pattern ends with trough or peak@
@assume ends with peak@


nrp=minc(nbt|nbp);

nrpp=nrp+1;

if nbt==nbp;
@ends with trough@
nrpp=nrp;
endif;

st[1:bcp5[1]]=ones(bcp5[1],1);
do while i<=nrpp;


st[bct5[i-1]+1:bcp5[i]]=ones(bcp5[i]-bct5[i-1],1);

i=i+1;
endo;

@test if no peaks equal no troughs@
@pattern then is pt,pt,pt etc and fill in from trough to end with ones@


if nbt==nrpp;

st[bct5[nrp]+1:nd]=ones(nd-bct5[nrp],1);
endif;

ret:

retp(st);
endp;


/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/


proc(1)= indicat(ddy);

local ff;
ff=0;
if ddy <= 0;
ff= 1;

endif;
retp(ff);
endp;


/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/


proc(0) = histo(x,bcp5,bct5,nd,nbp,nbt);     @ used to plot expansions @
local ntr, npk;
local td,pd,pdc,tdc,pde,tde,pda,tda, tdea, pdea;
local za,r,k,nr,nv,IR,r1,ag,in,x0,roe,yy;

if bcp5[1,1] .< bct5[1,1]; @ Peaks are first @

nr=nbt|nbp;
r=NBp-1;


k=1;
do while k<=r;
pdc=x[bct5[k,1]:bcp5[k+1,1],1];
?;

in = (pdc[rows(pdc),.] - pdc[1,.])/(rows(pdc)-1);
in = in*ones(rows(pdc),1);
x0=pdc[1,.];					       
roe=1;
yy=recserar(in,x0,roe);

pdc;
xy(seqa(1,1,rows(pdc)),pdc~yy);
wait;

pdc = pdc[2:rows(pdc),.] - pdc[1:rows(pdc)-1,.];
?;
@bar(0,pdc);
wait;@
k=k+1;
endo;

else;                      @ troughs are First@

r=NBt-1;
k=1;
do while k<=r;
pdc=x[bct5[k,1]:bcp5[k,1],1];
?;
in = (pdc[rows(pdc),.] - pdc[1,.])/(rows(pdc)-1);
in = in*ones(rows(pdc),1);
x0=pdc[1,.];					       
roe=1;
yy=recserar(in,x0,roe);

pdc;
xy(seqa(1,1,rows(pdc)),pdc~yy);
wait;

pdc = pdc[2:rows(pdc),.] - pdc[1:rows(pdc)-1,.];

?;
@bar(0,pdc);
wait;@
k=k+1;
endo;
endif;


endp;


/*************************************************************************************/
/*************************************************************************************/
/*************************************************************************************/

proc(4)=rawall(y);
local i,adf,adb,p1,jp,jt,ff,fis,fpt,j,ltp,jj;

@ finds initial turning points with restrictions in one step @

bcp5 = {};
bct5 = {};


adf=zeros(4,1);
i=turnphase+1;

do while i<=nd-turnphase;   @ to decide first  turning point @

   {fis,fpt} = isdate(y[i-turnphase:i+turnphase]);   
  
    if fis==1.and(fpt==1);
    bcp5= bcp5|i;
    goto proce;
    elseif fis==1.and(fpt==-1);   
    bct5=bct5|i;
    goto proce;
    endif;
        
    i=i+1;
endo;
proce:

if rows(bcp5)+rows(bct5) == 0;
nbp =0;nbt=0;bcp5=0;bct5=0;
goto nopt;
endif;


if fpt == 1;
j = bcp5; ltp = fpt;
i = j+1;
jj=0;
else;
j = bct5; ltp = fpt;
jj= bct5;
i = j+1;
endif;



do while i<=nd-turnphase;   @ find next turning point subject to restrictions @

   {fis,fpt} = isdate(y[i-turnphase:i+turnphase]);
  
  
    if rows(bcp5)+rows(bct5)<2;  @ needed as no full cycle - find second turning point @
    
        if fis==1.and(fpt==1).and(i-j>=phase).and(ltp==-1).and(indicat(y[i]-y[jj])==0); 
        bcp5= bcp5|i;
        ltp = fpt;
        j=i;
        elseif fis==1.and(fpt==-1).and(i-j>=phase).and(ltp==1).and(indicat(y[j]-y[i])==0),
        or(fis==1).and(fpt==-1).and(y[i]-y[j]<-thresh).and(ltp==1);
        bct5=bct5|i;
        ltp = fpt;
        jj=i;
        elseif fis==1.and(fpt==1).and(ltp==1).and(indicat(y[j]-y[i])==1); @ case of two peaks in a row @
        bcp5 = {};  @need to delete the last bcp5@
        bcp5= bcp5|i;
        ltp = fpt;
        j=i;
        elseif fis==1.and(fpt==-1).and(ltp==-1).and(indicat(y[i]-y[j])==1); @ case of two troughs in a row @
        bct5 = {};
        bct5=bct5|i;
        ltp = fpt;
        jj=i;
        endif;
        
    else;  @ all restrictions imposed - finds the rest of the turning points @
    
    
        if fis==1.and(fpt==1).and(i-jj>=phase).and(ltp==-1).and(i-j>=cycle).and(indicat(y[i]-y[jj])==0);
        bcp5= bcp5|i;
        ltp = fpt;
        j=i;
  
        goto next;
        
        elseif fis==1.and(fpt==1).and(ltp==-1).and(indicat(y[i]-y[jj])==0);  @ deals with situations were restrictions are breached  - peaks@
        
        
                if rows(bct5) > 1;
              
        
                        if (y[i] < y[j]).and(i-jj<phase ).or(y[i] < y[j]).and(i-jj<phase ).and(i-j<cycle)
                                .or(y[i] < y[j]).and(i-j<cycle)
                                        .or( (y[i] >= y[j]).and(y[jj]<=y[bct5[rows(bct5)-1]]).and(i-jj<phase ))
                                        .or( (y[i] >= y[j]).and(y[jj]<=y[bct5[rows(bct5)-1]]).and(i-jj<phase ).and(i-j<cycle));
                        
                                 @ effectively gets rid of this peak  @
                                 
                        elseif    (y[i] >= y[j]).and(y[jj]>y[bct5[rows(bct5)-1]]).and(i-jj<phase )
                                   .or((y[i] >= y[j]).and(y[jj]>y[bct5[rows(bct5)-1]]).and(i-j<cycle)) ;
                                     
                                   bct5 = trimr(bct5,0,1);
                                   
                                   if rows(bcp5) ==1;
                                    bcp5 = {};
                                   else;
                                    bcp5 = trimr(bcp5,0,1);
                                   endif;
                                                                                     
                                   bcp5= bcp5|i;
                                   ltp = fpt;
                                   j=i;
                                    
                        elseif (y[i] >= y[j]).and(y[jj]<=y[bct5[rows(bct5)-1]]).and(i-j<cycle);
                        
                                  if rows(bct5) == 2;
                                  bct5 = trimr(bct5,1,0) ; 
                                  else;
                                  bct5 = trimr(bct5,0,2)| bct5[rows(bct5)] ; 
                                  endif; 
                                  
                                  if rows(bcp5) == 1;
                                    bcp5 = {};
                                  else;
                                    bcp5 = trimr(bcp5,0,1);
                                  endif;                                    
                                      
                                  bcp5= bcp5|i;
                                  ltp = fpt;
                                  j=i;
                        
                         endif;
                         
                elseif rows(bct5) == 1;
                
                
                        if (y[i] < y[j]).and(i-jj<phase ).or(y[i] < y[j]).and(i-jj<phase ).and(i-j<cycle)
                                .or(y[i] < y[j]).and(i-j<cycle).or( (y[i] >= y[j]).and(i-jj<phase ))
                                        .or( (y[i] >= y[j]).and(i-jj<phase ).and(i-j<cycle));
                        
                                 @ effectively gets rid of this peak @
                                
                                              
                        elseif (y[i] >= y[j]).and(i-j<cycle);
                        
                                  bcp5 = {};
                                      
                                  bcp5= bcp5|i;
                                  ltp = fpt;
                                  j=i; 
                        
                         endif;
                                         
                         
                endif;                   
          
            
                  
        elseif fis==1.and(fpt==-1).and(i-j>=phase).and(ltp==1).and(i-jj>=cycle).and(indicat(y[j]-y[i])==0),
        or(fis==1).and(fpt==-1).and(y[i]-y[j]<-thresh).and(ltp==1);
        bct5=bct5|i;
        ltp = fpt;
        jj=i; 
        
        goto next;
        
        elseif fis==1.and(fpt==-1).and(ltp==1);  @ deals with situations were restrictions are breached  - troughs @
        
                
            if rows(bcp5) > 1;    
                
                
                 if (y[i] > y[jj]).and(i-j<phase ).or(y[i] > y[jj]).and(i-j<phase ).and(i-jj<cycle)
                        .or(y[i] > y[jj]).and(i-jj<cycle)
                        .or( (y[i] <= y[jj]).and(y[j]>=y[bcp5[rows(bcp5)-1]]).and(i-j<phase ))
                        .or( (y[i] <= y[j]).and(y[j]>=y[bcp5[rows(bcp5)-1]]).and(i-j<phase ).and(i-jj<cycle));  
              
                       @ effectively gets rid of this trough @
                       
                       
                elseif (y[i] <= y[jj]).and(y[j]<y[bcp5[rows(bcp5)-1]]).and(i-j<phase )
                       .or((y[i] <= y[j]).and(y[j]<y[bcp5[rows(bcp5)-1]]).and(i-jj<cycle)) ;
                       
                         
                       bcp5 = trimr(bcp5,0,1);
                       
                       if rows(bct5) ==1;
                        bct5 = {};
                       else;
                        bct5 = trimr(bct5,0,1);
                       endif; 
                       
                       bct5= bct5|i;
                       ltp = fpt;
                       jj=i;
                
                elseif(y[i] <= y[jj]).and(y[j]>=y[bcp5[rows(bcp5)-1]]).and(i-jj<cycle);
                   
                      if rows(bcp5) == 2;
                        bcp5 = trimr(bcp5,1,0) ; 
                      else;
                         bcp5 = trimr(bcp5,0,2)| bcp5[rows(bcp5)] ; 
                      endif; 
                      
                      if rows(bct5) == 1;
                        bct5 = {};
                      else;
                        bct5 = trimr(bct5,0,1);
                      endif;  
                      
                      bct5= bct5|i;
                      ltp = fpt;
                      jj=i;
                
                
                endif;
                
            
            elseif rows(bcp5) == 1;
            
                 if (y[i] > y[jj]).and(i-j<phase ).or(y[i] > y[jj]).and(i-j<phase ).and(i-jj<cycle)
                        .or(y[i] > y[jj]).and(i-jj<cycle)
                        .or( (y[i] <= y[jj]).and(i-j<phase ))
                        .or( (y[i] <= y[j])and(i-j<phase ).and(i-jj<cycle));  
              
                       @ effectively gets rid of this trough @
                
                 elseif(y[i] <= y[jj]).and(i-jj<cycle);
                   
                      bct5 = {};
               
                      bct5= bct5|i;
                      ltp = fpt;
                      jj=i;
                endif;                      
                
            endif;    
        
           
        
        
        
        
        elseif fis==1.and(fpt==1).and(ltp==1).and(indicat(y[j]-y[i])==1);       @ case of two peaks in a row @;
            if rows(bcp5) == 1;
            bcp5 = {};
            else;
            bcp5 = trimr(bcp5,0,1);  @need to delete the last bcp5@
            endif;
        bcp5= bcp5|i;
        ltp = fpt;
        j=i;
        
        elseif fis==1.and(fpt==-1).and(ltp==-1).and(indicat(y[i]-y[jj])==1);    @ case of two troughs in a row @;
         
            if rows(bct5) == 1;
            bct5 = {};
            else;
            bct5 = trimr(bct5,0,1);  @need to delete the last bcp5@
            endif;
        bct5=bct5|i;
        ltp = fpt;
        jj=i;
        endif;
        
    endif;
    
    
    next:
   @ i~bct5[rows(bct5),.];@
    i=i+1;
endo;

if rows(bcp5)+rows(bct5) == 1;
nbp =0;nbt=0;bcp5=0;bct5=0;
goto nopt;
endif;



if (bcp5[1,.] > bct5[1,.]);
 
    j = bct5[1,.]-1;
    do while j > 1;
        if y[j,.]<y[bct5[1,.],.];
            if rows(bct5) ==1;
            bct5 = {};
            goto finb;
            else;       
            bct5 = trimr(bct5,1,0);
            goto finb;  
            endif;
        endif;
        
        j = j-1;
    endo;

endif;

if (bct5[1,.] > bcp5[1,.]);
 
    j = bcp5[1,.]-1;
    do while j > 1;
        if y[j,.]>y[bcp5[1,.],.];
            if rows(bcp5) ==1;
            bcp5 = {};
            goto finb;
            else;     
            bcp5 = trimr(bcp5,1,0);
            goto finb;
            endif;  
        endif;
        
        j = j-1;
    endo;

endif;

finb:

if rows(bcp5)+rows(bct5) == 1;
nbp =0;nbt=0;bcp5=0;bct5=0;
goto nopt;
endif;


nbp = rows(bcp5);
nbt = rows(bct5);

if (bcp5[nbp,.] > bct5[nbt,.]);
 
    j = bcp5[nbp,.]+1;
    do while j < nd;
        if y[j,.]>y[bcp5[nbp,.],.];
            if rows(bcp5)>2;
            bcp5 = trimr(bcp5,0,1);
            else;
            bcp5 = {};
            endif;
        goto fin;
        endif;
        
        j = j+1;
    endo;

endif;


if (bcp5[nbp,.] < bct5[nbt,.]);
    
    j = bct5[nbt,.]+1;
    do while j < nd;
        if y[j,.]<y[bct5[nbt,.],.];
            if rows(bct5)>2;
            bct5 = trimr(bct5,0,1);
            else;
            bct5 = {};
            endif;
   
        goto fin;
        endif;
        
        j = j+1;
    endo;

endif;







fin:
nbp = rows(bcp5);
nbt = rows(bct5);
nopt:
retp(bcp5,bct5,nbp,nbt);
endp;




proc(2)= isdate(y);

@ is this point a turning point and is it peak or trough @
@ use turnphase to determine if peak or trough @

local ff,adf,p1,i,fpt,fis,q,qq;
adf=zeros(turnphase*2,1);

i = turnphase+1;

fis = 0;
fpt = 0;

    q = 1;
    qq = 1;
    
    do while q <= turnphase*2+1;
    
        if q/=i;
        adf[qq] = indicat(y[q]-y[i]); 
        qq = qq+1;
        
        endif;
    
    q = q+1;
    endo;
     
     
    p1 = sumc(adf);

    
    if p1==turnphase*2;
    fis = 1;
    fpt = 1;
  
    endif;
    
    
    
    q = 1;
    qq = 1;
    
    do while q <= turnphase*2+1;
        if q/=i;
        adf[qq] = indicat(y[i]-y[q]); 
        qq = qq+1;
        
        endif;
    
    q = q+1;
    endo;
          
    p1 = sumc(adf);
    
    if p1==turnphase*2;
    fis = 1;
    fpt = -1;
    
    endif;


retp(fis,fpt);
endp;



/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/

/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/


@ BELOW HERE IS MODEL SIMULATION CODE    @







































/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/
/******************************************************************************************************/

/* include model simulation code here */



proc(1) = usar; 


local mn,st,init,e,u,ht,roe,x0,xx,i,xxm;

@Simulating from an AR(3)@

@will drop first 500 observations so computes 500 more than needed@
nd=nd+500;

xxm=zeros(nd,1);
xx=zeros(nd,1);
e=rndn(nd,1);

st=.01;
u=e*st;

i = 4;
xxm[1] = 0;
xxm[2] = 0;
xxm[3]=0;

do while i <= nd;

xxm[i] = 1.55*xxm[i-1] - .675*xxm[i-2]+  u[i];
@if random walk remove @ 
@xx[i]=xx[i-1]+xxm[i-1];@
xx[i]=xxm[i];
i = 1+i;
endo;


x=xx[501:nd];
nd=nd-500;
retp(x);
endp;