/* An example for calling logpower SAS macro
%include 'logpower.sas';
%logpower(n=20, mu=0.5, sigma=0.8, M0=1.1, alpha=0.05);

run;

/**********************************************************************
 NAME: logpower.sas

 This program computes the power of the generalized test for testing
 the mean of a lognormal distribution with parameter mu and sigma.
 mu = population mean of the logged data and  sigma = population
 standard deviation of the logged data. For a given sample size n, mu
 and sigma, this program computes the power for testing

 H_0: M >= M0    vs    H_a: M < M0      (*)

 where M = mu + sigma^2/2 and M0 is a specified value of M.
 Note that exp(M) is the mean of the lognormal distribution. Testing
 (*) is equivalent to testing

 H_0: exp(M) >= exp(M0)    vs    H_a: exp(M) < exp(M0)

  INPUT:
            n = sample size
           mu = population mean of the logged data
        sigma = population standard deviation of the logged data
           M0 = a specified value of M
        alpha = test level
           m1 = number of pwoers for generalized test. Default = 2500
           m2 = number of pavalues: Default = 5000

 OUTPUT: power
************************************************************************/


/* program of power of the test of lognormal mean */

%macro logpower(n=, mu=, sigma=, M0=, alpha=, m1=2500, m2=5000);

options ls = 64 ps = 45 nodate nonumber;
title 'Output of power of test of lognormal mean';
proc iml;

n=&n;
ul=&mu;
sigmal=&sigma;
mean0=&M0;
alpha=&alpha;
m1 = &m1;
m2 = &m2;

reset name=noname center=nocenter;

df = n-1.0;
seed1 = int(time());
seed2 = int(time()+12345);
seed3 = int(time()+13456);
seed4 = int(time()+19867);
power = 0.0;
do i= 1 to m1;
   z1 = rannor(seed1);
   v1 = 2.0*rangam(seed2,df/2.0);
   var = sigmal**2*v1/df;
   s = sqrt(var);
   xbar = ul +z1*sigmal/sqrt(1.0*n);
   const1 = s*sqrt(df)/sqrt(n);
   const2 = var*df*0.5;

   pvalue = 0.0;

   do j = 1 to m2;
      z2 = rannor(seed3);
      v2 = 2.0*rangam(seed4,df/2.0);
      stat = xbar + z2*const1/sqrt(v2)+const2/v2;
      if stat > mean0 then  pvalue = pvalue + 1.0;
   end;
   pvalue = pvalue/m2;

   if pvalue < alpha then  power = power + 1.0;
end;
power = power/m1;

print 'n=' n[format=2.0] ',' 'ul=' s[format=6.3] ','
   'sigmal=' sigma[format=6.3] ',' 'mean0=' mean0[format=6.3];
print 'The power for right-tail test is: ' power;

quit;

%mend logpower;