goptions reset=all; OPTIONS NODATE NONUMBER linesize=100; DM 'CLEAR LOG'; DM 'CLEAR OUTPUT';

%Let Loc = C:\SAS Files\Monte Carlo Simulations ; Libname Loc "&Loc\Dissertation\Phase I\";

/* ------------------------------------------------------------------------ --------------- */ %Let L3 = .007 ; /*--- INPUT SEED FOR LAMBDA 3 ---*/ %Let L4 = &L3 ; /*--- INPUT SEED FOR LAMBDA 4 ---*/ %Let Mean = 0; %Let Var = 1; %Let Inc = .5; %Let minSkew = -2; %Let minKurt = -1; %Let maxSkew = 2; %Let maxKurt = 5.5; %Let CV =.00000001; /* DEFAULT CONVERGENCE CRITERION */ %Let MAXITER=500; /* DEFAULT MAXIMUM ITERATIONS */ /* ------------------------------------------------------------------------ --------------- */ Title 'Fitting Distributions with the GLD via the Method of Moments' ; %Macro FindLambda(a1=, a2=, a3=, a4=); PROC IML; START Moments ; MEAN = &a1; VAR = &a2; SKEW = &a3; %Let S = &a3 ; KURT = &a4; IF SKEW < 0 THEN DO; SKEW = abs(SKEW) ; %Let S = %sysfunc(abs(%sysevalf(&a3))) ; END; FINISH;

START Valid; IF %sysfunc(min(&a4,%sysevalf(&S**2 - 2)))=&a4 THEN DO; minvalue= %sysevalf(&S**2 - 2); PRINT "INVALID REGION: Skew=" &S " Kurt=" &a4; PRINT "You must increase Kurtosis to more than: " minvalue; ABORT; END; FINISH;

START NEWTON; /* -------- ---------------------------- */ RUN FUN; /* EVALUATE FUNCTION AT STARTING VALUES */ DO ITER=1 TO &MAXITER /* ITERATE UNTIL MAXITER ITERATIONS */ WHILE(MAX(ABS(F))>&CV); /* OR CONVERGENCE */ RUN DERIV; /* EVALUATE DERIVATIVES IN J */ DELTA=-SOLVE(J,F); /* SOLVE FOR CORRECTION VECTOR */ Lambda=Lambda+DELTA; /* THE NEW APPROXIMATION */ RUN FUN; /* EVALUATE THE FUNCTION */ END; /* ------------------------------------ */ FINISH;

/*---USER-SUPPLIED FUNCTION EVALUATION---*/ START FUN; L3=Lambda[1]; L4=Lambda[2]; /* EXTRACT THE VALUES */

A = 1/(1+L3) - 1/(1+L4); B = 1/(1+2*L3) + 1/(1+2*L4) - 2*BETA(1+L3,1+L4); C = 1/(1+3*L3) - 1/(1+3*L4) - 3*BETA(1+2*L3,1+L4) + 3*BETA (1+L3,1+2*L4); D = 1/(1+4*L3) + 1/(1+4*L4) - 4*BETA(1+3*L3,1+L4) + 6*BETA (1+2*L3,1+2*L4) - 4*BETA(1+L3,1+3*L4);

L2 = SQRT((B-A**2)/&a2); L1 = &a1-A/L2;

/* The first part of F estimates SKEWNESS the second part estimates KURTOSIS */ F= (((B-A**2)**(-1.5))*(C-3*A*B+2*A**3) - SKEW)//

(((B-A**2)**(-2))*(D-4*A*C+6*(A**2)*B-3*A**4) - (KURT+3));/* 3 is added to */ FINISH; /* normalize kurtosis*/

/*---user-supplied derivatives of the function---*/ START DERIV; dA_L3 = -(1+L3)**(-2); dB_L3 = -2*(1+2*L3)**(-2) - 2*BETA(1+L3,1+L4)*(DIGAMMA (1+L3)-DIGAMMA(2+L3+L4)); dC_L3 = -3*(1+3*L3)**(-2) - 6*BETA(1+2*L3,1+L4)*(DIGAMMA (1+2*L3)-DIGAMMA(2+2*L3+L4)) +3*BETA(1+L3,1+2*L4)*(DIGAMMA(1+L3)-DIGAMMA (2+L3+2*L4)); dD_L3 = -4*(1+4*L3)**(-2) - 12*BETA(1+3*L3,1+L4)*(DIGAMMA (1+3*L3)-DIGAMMA(2+3*L3+L4)) +12*BETA(1+2*L3,1+2*L4)*(DIGAMMA(1+2*L3)-DIGAMMA (2+2*L3+2*L4)) -4*BETA(1+L3,1+3*L4)*(DIGAMMA(1+L3)-DIGAMMA (2+L3+3*L4));

dA_L4 = (1+L4)**(-2); dB_L4 = -2*(1+2*L4)**(-2) - 2*BETA(1+L3,1+L4)*(DIGAMMA (1+L4)-DIGAMMA(2+L3+L4)); dC_L4 = 3*(1+3*L4)**(-2) - 3*BETA(1+2*L3,1+L4)*(DIGAMMA (1+L4)-DIGAMMA(2+2*L3+L4)) +6*BETA(1+L3,1+2*L4)*(DIGAMMA(1+2*L4)-DIGAMMA (2+L3+2*L4)); dD_L4 = -4*(1+4*L4)**(-2) - 4*BETA(1+3*L3,1+L4)*(DIGAMMA (1+L4)-DIGAMMA(2+3*L3+L4)) +12*BETA(1+2*L3,1+2*L4)*(DIGAMMA(1+2*L4)-DIGAMMA (2+2*L3+2*L4)) -12*BETA(1+L3,1+3*L4)*(DIGAMMA(1+3*L4)-DIGAMMA (2+L3+3*L4));

da3L3 = (B-A**2)**(-1.5)*(dC_L3-3*(A*dB_L3+B*dA_L3) + 6* (A**2)*dA_L3) + (C-3*A*B+2*A**3)*(-1.5*(B-A**2)** (-2.5)*(dB_L3-2*A*dA_L3)) ; da3L4 = (B-A**2)**(-1.5)*(dC_L4-3*(A*dB_L4+B*dA_L4) + 6* (A**2)*dA_L4) + (C-3*A*B+2*A**3)*(-1.5*(B-A**2)** (-2.5)*(dB_L4-2*A*dA_L4)) ;

da4L3 = (B-A**2)**(-2)*(dD_L3-4*(A*dC_L3+C*dA_L3) + 6* (A**2*dB_L3+2*A*B*dA_L3) - 12*(A**3)*dA_L3) + (D-4*A*C+6*(A**2)*B-3*A**4)*(-2* (B-A**2)**(-3)*(dB_L3-2*A*dA_L3)) ; da4L4 = (B-A**2)**(-2)*(dD_L4-4*(A*dC_L4+C*dA_L4) + 6* (A**2*dB_L4+2*A*B*dA_L4) - 12*(A**3)*dA_L4) + (D-4*A*C+6*(A**2)*B-3*A**4)*(-2* (B-A**2)**(-3)*(dB_L4-2*A*dA_L4)) ;

J=((da3L3)||(da3L4))// ((da4L3)||(da4L4)); /* EVALUATE JACOBIAN */ FINISH;

START NEGATIVE; If L2>0 & L3<0 & L4<0 then Do; L2 =-L2; T = L3; L3 = L4; L4 = T; end; If &a3 < 0 then DO; SKEW = -1*SKEW; T = L3; L3 = L4; L4 = T; L1 = L1 + 2*A/L2; end; FINISH;

DO; Lambda={&L3, &L4}; /* starting values */ RUN Moments ; RUN Valid ; RUN NEWTON; RUN NEGATIVE;

PRINT "MOMENTS" ; PRINT MEAN VAR SKEW KURT ; PRINT "Estimates of Lambdas 1 to 4." ; PRINT L1 L2 L3 L4 ;

CREATE Work.Lambdas var{MEAN VAR SKEW KURT L1 L2 L3 L4}; APPEND; CLOSE Work.Lambdas ; END ; run; Title; %Mend FindLambda;

%MACRO AppendData; *%put Does file exist (1=Yes, 0=NO)? %sysfunc(exist(work.Lambdas)); %if %sysfunc(exist(work.Lambdas))=1 %then %do ; /* Appends results of each run to _Data. */ Proc Append Base=Loc.Lambdas Data = Work.Lambdas ; run ;

proc datasets ; delete Lambdas; run; %end; %else %put Temporary file Lambdas does not exist in the Work folder.; %Mend;

%MACRO Lambdas(Kurt=); %DO %while(%sysevalf(&minSkew<=&maxSkew)); %Let Kurt = &minKurt; %DO %while(%sysevalf(&Kurt<=&maxKurt)); %put skew=%quote(&minSkew) kurt=%quote(&Kurt); %FindLambda(a1=&Mean, a2=&Var, a3=&minSkew, a4=&Kurt); %AppendData; %let Kurt=%sysevalf(&Kurt+&Inc); %END; %let minSkew=%sysevalf(&minSkew+&Inc); %END; %MEND; %Lambdas(Kurt=&minKurt); DATA Loc.Lambdas; Retain _Dist_ Distribution N Mean Var Skew Kurt sd_Skew sd_Kurt L1 L2 L3 L4; Format skew kurt 4.1; set Loc.Lambdas; _Dist_ = _N_; run; QUIT;