/*
The macro %DIFFRR constructs confidence intervals around the absolute and relative effect estimates from 
segmented linear regressions using Delta method. We run the selected final segmented autoregressive error model, 
routing the output to a pre-designated file and writing to data sets the parameter estimates and covariance matrix 
obtained through maximum likelihood estimation, assuming as given the autocorrelations estimated by the procedure.  
Then we choose the maximum likelihood estimates, if available.  Otherwise, we use OLS estimates if maximum likelihood
estimation does not converge. One can use Yule-Walker or other types of estimates by changing the corresponding codes.
Finally, we use the IML procedure to calculate the expectation and variance of the absolute difference as well as 
the relative change using the multivariate delta method for the latter.  
Normal distribution (based on large sample theory) is assumed to construct the final confidence intervals.
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 seven required parameters for the macro:

In - the name of the time series data set to be used to run the final model
Out - the output data set with the information created by the macro
Outcome - the dependent variable of the time series regression model
Predictors - the independent variables in the final time series regression model 
Lag - the lag term for the final time series regression model
Variablewithout - the values of the variables in the final model at the time point of interest, without the intervention
Variablewith - the values of the variables in the final model at the time point of interest, with the intervention
*/
%macro DIFFRR(in,out,outcome,predictors,lag,coefwithout,coefwith,siglevel);
%macro sort(dsn,ods,by,where,keep);

 proc sort data=&dsn(keep=&keep) out=&ods;
      by &by;
      where &where;
 run;

%mend sort;
filename routed 'c:\zzzz.txt';
proc printto print=routed new;
run;
/*we output the parameter estimates and covariate matrix of the final ML model that assumes the autregressive parameters given*/
proc autoreg data=&in outest=est covout;
 ods output Autoreg.Model1.FinalModel.CovB=work.afinalcov∈
 ods output Autoreg.Model1.FinalModel.ParameterEstimatesGivenAR=work.finalparam&in(keep=estimate variable);
 ods output Autoreg.Model1.OLSEst.CovB=work.olscov&in(drop=model);
 ods output Autoreg.Model1.OLSEst.ParameterEstimates=work.olsparam&in(keep=estimate variable);
 model &outcome=&predictors/
 method=ml nlag=&lag backstep dwprob covb;
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 ne 1 %then %do;
    data afinalcov∈set afinalcov&in(keep=rowname intercept &predictors);
     j=_n_;
    run;
    %sort(afinalcov&in,tfinalcov&in,rowname j);
    data finalcov&in(where=(substr(rowname,1,2) ne 'AR'));set tfinalcov∈by rowname;if last.rowname;run;
    %sort(finalcov&in,finalcov&in,j);

    data d1;set finalcov&in(drop=j);run;
    data d2;set finalparam∈run;
%end;
%else %do;data d1;set olscov∈run;
          data d2;set olsparam∈run;
%end;


/*we use the output of the parameter estimates of the final model*/
proc transpose data=d2 out=param(drop=_name_);
     var estimate;
	 id variable;
run;

proc iml;
      use d1;
      read all into x;
	  use param;
      read all into est; 
	  c={&coefwithout};
	  d={&coefwith};
	  dis=d-c;
	  ywithout=c*est`;
	  ywith=d*est`;
	  varwithout=c*x*c`;
	  varwith_dis=dis*x*dis`;
	  cov=c*x*dis`;
	  absdiff=(ywith-ywithout);  
	  relchange=(ywith-ywithout)/ywithout;  
	  y=j(1,10);
	  varabsdiff=dis*x*dis`;
	  varrelchange=(1/(ywithout**2))*varwith_dis+(absdiff**2/(ywithout**4))*varwithout-
                    2*(absdiff/(ywithout**3))*cov;
	  hfcit=1-(&siglevel)/2;
	  crt=quantile('NORMAL',hfcit);
      lowerabsdiff=absdiff-crt*sqrt(varabsdiff);
	  upperabsdiff=absdiff+crt*sqrt(varabsdiff);
      lowerrelchange=relchange-crt*sqrt(varrelchange);
	  upperrelchange=relchange+crt*sqrt(varrelchange);
	  y[1,1]=ywithout;y[1,2]=ywith;
      y[1,3]=absdiff;y[1,4]=varabsdiff;y[1,5]=lowerabsdiff;y[1,6]=upperabsdiff;
      y[1,7]=relchange;y[1,8]=varrelchange;y[1,9]=lowerrelchange;y[1,10]=upperrelchange;
      create variance from y;
	  append from y;
quit;
data &out;set variance;rename 
   col1=y1 col2=y2 col3=absdiff col4=varabsdiff col5=lowerabsdiff col6=upperabsdiff
   col7=relchange col8=varrelchange col9=lowerrelchange col10=upperrelchange;
run;
%mend DIFFRR;