%macro trap_area_Z(LOW=-1.645, HIGH=1.645, npts_lo=10, npts_hi=10, npts_by=2,
   fout=C:\Users\baileraj\Desktop\est3.out,
   print_est=FALSE, print_pts=FALSE, display_graph=FALSE, ODS_on=FALSE);

/*  ======================================================================== 
Purpose: Estimate P{LOW<Z< HIGH) using the trapezoidal rule
Macro variables:
LOW, HIGH: interval of interest
NPTS_LO, NPTS_HI, NPTS_BY: # function values evaluated in area calc.
FOUT: output data file containing area estimate for each NPTS value
PRINT_EST: print area estimate
PRINT_PTS: print points/nodes {x1-xn} + function values {f(x1)-f(xn)}   
DISPLAY_GRAPH: generate PROC GPLOT with function values
ODS_ON: generate ODS RTF output
======================================================================== 
*/

%do npts = &npts_lo %to &npts_hi %by &npts_by; *loop over npts values;

data trapper;
  file "&fout" MOD;
  trapsum = 0;
  array x_value(&npts) x1-x&npts;
  array f_value(&npts) y1-y&npts;

  low = &LOW;
  high = &HIGH;
  incr = (high-low)/( &npts -1);
  pi = arcos(-1);
  do i= 1 to &npts;
    x_value[i] = low + incr*(i-1);
    f_value[i] = (1/sqrt(2*pi))*exp(-x_value[i]*x_value[i]/2);

    if i=1 or i=&npts then trapsum = trapsum + f_value[i]/2;
	else trapsum = trapsum + f_value[i];
  end;
  area_est = trapsum*incr;
  output;
 
put "est. P(&LOW < Z < &HIGH) =" area_est "(based on &NPTS points)";
run;

%if %upcase(&ODS_ON)=TRUE %then %do;
   ods rtf bodytitle file="ods_out.rtf";
%end;

%if %upcase(&print_est)=TRUE %then %do;
proc print data=trapper;
  title "Trapezoidal Rule Area Estimate for P(&LOW<Z<&HIGH)";
  title2 "(based on &NPTS equally spaced points)";
  var low high incr area_est;
run;
%end;

data trapper2;
   set trapper;
   array x_value(&npts) x1-x&npts;
   array f_value(&npts) y1-y&npts;
   do ii=1 to &npts;
     xout = x_value[ii];
     yout = f_value[ii];
     output;
   end;
run;

%if %upcase(&print_pts)=TRUE %then %do;
proc print data=trapper2;
title "Interpolation Points for Trapezoidal Rule";
   var ii low high incr area_est xout yout;
run;
%end;

%if %upcase(&display_graph)=TRUE %then %do;
ods graphics on;
proc sgplot data=trapper2;
title "Plot of function values vs. x-values";
  scatter x=xout y=yout;
run;
ods graphics off;
%end;

%if %upcase(&ODS_ON)=TRUE %then %do;
ods rtf close;
%end;

%end;   * of loop over npts values;
%mend trap_area_Z;

options mprint;
%trap_area_Z(LOW=0,HIGH=1.96,npts_lo=5,npts_hi=25,
npts_by=5, display_graph=TRUE)