/*
The macro %BOOTS constructs confidence intervals around the absolute and relative effect estimates from 
segmented linear regressions using bootstrapping method. We simulate a time series based on the final selected model 
at each iteration. We route the output to a pre-designated file so that we can exclude iterations where the maximum 
likelihood method does not converge. We rerun the segmented time series regression to obtain estimates for each valid
iteration. We summarize the distribution of all estimates from the simulations and construct the confidence interval
from data set simulation results with flag delete=0.

1.	Identify a temporary file to which the SAS output will be routed by changing the path in the first statement.  
(For example: filename routed 'c:\zzzz.txt';)

2.	Specify the following required parameters for the macro:

Baseint - the estimated baseline intercept from the final model
Baseslope - the estimated baseline slope from the final model
Interaint - the estimated level change comparing the intervention period to the baseline from the final model
Interaslope - the estimated trend change comparing the intervention period to the baseline from the final model 
Nperiod - the number of total time points
Sigmas - the estimated variance from the final model
Ar - the estimated autocorrelation terms from the final model 
Iter - the number of simulation iterations 
Intma - the value of the starting time point of the first policy
Modar - the maximum allowable autocorrelations will be used in each iteration of segmented regression desired by the investigators
Fieldb - the list of significant terms from the final model
*/
data logall;intercept=0;time=0;
     intervention1=0;timeafterintervention1=0;
     delete=1;
run;
%macro BOOTS(baseint,baseslope,interaint,interaslope,
   nperiod,sigmas,ar,iter,intma,modar,fieldb);
   %macro numa;%do i=1 %to &iter;logd&i%end;%mend numa;
   data xbeta(drop=beta0-beta3);
        beta0=&baseint;beta1=&baseslope;beta2=&interaint;beta3=&interaslope;
		do time=1 to &nperiod;
	       intervention1=(time >=&intma);
	       timeafterintervention1=max(0,time-&intma+1);
           mu2=beta0+beta1*time+beta2*intervention1+beta3*timeafterintervention1;
		   output;
	    end;
   run;
   %do k=1 %to &iter;
   proc iml; 
        reset noname; 
        use xbeta var {mu2};
        read all into z2;
        start armasim(y,n,phi,theta,seed,sigma_2); 
           p = ncol(phi)-1;q = ncol(theta)-1; 
           y = normal(j(1,n+q,seed));y = sqrt(sigma_2)*y;
        /*----  Pure MA or white noise  ----*/ 
        if p=0 then y=product(theta,y)[,(q+1):(n+q)]; 
        else do;   /*----  Pure AR or ARMA  ----*/ 
         /*----  Get the autocovariance function  ----*/ 
              call armacov(gamma,cov,ma,phi,theta,p); 
              if gamma[1]<0 then do; 
                 print 'ARMA parameters not stable.'; 
                 print 'Execution terminating.'; 
                 stop; 
              end; 
         /*----  Form covariance matrix  ----*/ 
              gamma=toeplitz(gamma); 
         /*----  Generate covariance between initial y and  ----*/ 
         /*----  initial innovations                        ----*/ 
              if q>0 then do; 
                 psi = ratio(phi,theta,q); 
                 psi = hankel( psi[,-( (-q):(-1) ) ] ); 
                 m = max(1,(q-p+1)); 
                 psi = psi[ -((-q):(-m)) ,]; 
                 if p>q then psi = j(p-q,q,0)//psi; 
                 gamma = ( gamma||psi ) // ( psi`||i(q) ); 
              end; 
         /*----  Use Cholesky root to get startup values  ----*/ 
              gamma = root(gamma); 
              startup = y[,1:(p+q)]*gamma; 
              e = y[,(p+q+1):(n+q)]; 
         /*----  Generate MA part  ----*/ 
              if q>0 then do; 
                 e = startup[,(p+1):(p+q)]||e; 
                 e = product(theta,e)[,(q+1):(n+q-p)]; 
              end; 
              y = startup[,1:p]; 
              phi1 = phi[, -(-(p+1):(-2)) ]`; 
         /*----  Use difference equation to generate  ----*/ 
         /*----  remaining values                     ----*/ 
              do ii=1 to n-p; 
                 y = y||( e[,ii]-y[,ii:(ii+p-1)]*phi1 ); 
              end; 
          end; 
              y = y`; 
       finish;     /*----  ARMASIM  ----*/ 
 
       seed=%eval(&k*&iter);
	   sample=%eval(&nperiod+200);
	   sam=sample+1-%eval(&nperiod);
	   run armasim(y,sample,{&ar},{1},seed,&sigmas);	   
       ya=y[sam:sample];
       yf=ya+z2;
       create ts from yf;
       append from yf;
  quit; 

  data dd;
       merge xbeta ts;
       /*if &phaseup ne . then do;if &phaseup>=time>=&phaselow then col1=.;end;
  	   else do;end;*/
  run;

  filename routed 'c:\zzzz.txt';
  proc printto print=routed new;
  run;
  proc autoreg data=dd outest=logd(keep=intercept &fieldb);
	   model col1=&fieldb/backstep method=ml nlag=&modar;	
  run;
 
  proc printto;run;

  data test;
     infile routed pad end=eof;
     input word $char256.;
     retain ml mle 0;
	 if index(word,'Maximum Likelihood Estimates') > 0 then ml=1;
	 if index(word,'ERROR: The estimation did not converge') >0 then
	    mle=1;
     if eof then call symput('ml',compress(ml)) ;
     if eof then call symput('mle',compress(mle)) ; 
  run; 

  %if &ml=1 and &mle=1 %then %do;
      data logd;intercept=0;time=0;
           intervention1=0;timeafterintervention1=0;
           delete=1;
      run;
	  proc append base=logall data=logd;run;
  %end;
  %else %if &ml=1 and &mle ne 1 %then %do;
      data logd;set logd;delete=0;
      proc append base=logall data=logd;run;
  %end;
  %else %do;
      data logd;set logd;delete=0;
      proc append base=logall data=logd;run;
  %end;
 
  %end;
  data simulationresults;set logall;run;
  
%mend BOOTS;