/* INTRO: This program calculates critical points for a multivariate normal with general variance covariance structure. These critical values can be used to perform MCC nd MCB inference per Horrace and Schmidt (1997). The critical values are for two-sided confidence interval construction, but with a minor modification (See below), can be adapted to one-sided values. For clarification of any aspect of this program contact Bill Horrace at whorrace@u.arizona.edu. OVERVIEW: The program generates a large number of k-dimensional multivariate normal random variables; finds the absolute value of the max of the k dimensions; ranks them; then find the (1-alpha)%ile of the order statistic. Rather that creating and ranking a vector of dimension equal to the simulation sample size, this program only generates the alpha*(simulation sample size) upper portion of the ranked vector. That is, for a simulation sample size of 1,000,000 and alpha = 0.05, this program creates and ranks a vector of size 0.05*1,000,000 = 50,000, rather than a vector of size 1,000,000. This saves alot of time when the sample size is large. To use this program the dos subdirectory that contains this program must contain a variance matrix for the parameters of interest before differencing. For instance if g' = [g(1) g(2) ... g(k)] are the parameters that will be differenced with the "best", [max g(i)-g(1), max g(i)-g(2),...max g(i)-g(k)], then the variance matrix is var(g). This matrix must be of dimension (kxk) and must be a gauss matrix called VARALPH.FMT. The user controls several variables that are defined in the initialization (INIT:) subroutine.... K is dimension of the multivariate normal distribution. For MCB it is the number of simultaneous comparisons [max g(i)-g(1),...max g(i)-g(k)]. ALPHA = (1 - confidence level). e.g. alpha = 0.05 for 95% confidence intervals. CTRL = the index of the control. MCB with general covariance structure is just k MCC's, so to perform MCB you will need to run this program k-times to generate k critical values. Alternatively you can modify this program to loop through the variable "ctrl" = 1, ..., k. SAMPLE = the number of times the multivariate normal random variable is sampled to generate the critical value. High SAMPLE means greater accuracy of the simulation but slower run time. TIMES = the number of multivariate normal variables that are pulled in each iteration of the simulation. That is, instead of pulling 1 realization of a multivariate normal SAMPLE times. It may be more economical to pull TIMES realizations of a multivariate SAMPLE/TIMES times. Obviously, SAMPLE/TIMES should be an integer. You may want to experiment with what will be the fastest TIMES for a given SAMPLE and k. HINTS: You may want to experiment with the TIMES value to see what works best for a given value of k. I have found that for k < 160, TIMES = 4000 is pretty economical. For 160 < k < 400, TIMES = 2000 is good. When experimenting with this it is a good idea to add a print statement in the SAMPLEUP loop. This way Gauss will send output to the screen each time it pulls TIMES realizations of the random variables. You can time how long it takes to get through a portion of the simulation sample size for different values of TIMES and optimize the program. Also, the print statements will allow you to stop the program with a "control-break" at the end of each experiment if you are running in DOS. Once you have decided on the optimal TIMES, remove the print statement to optimize performance. ACCURACY: the critical values are themselves random variables (unfortunately), so they have some associated sampling variability. I have found that with SAMPLE = 1,000,000 the critical values generated are "accurate" to 2 decimal places. I haven't experimented with other values of SAMPLE. ONE SIDED CRITICAL VALUES: This program generates two sided critical values by taking the absolute values of the realization of the generated random variables. To generate one-sided critical values remove the "abs(.)" from the statement "maxabn = maxc(abs(n))" in the "generate" subroutine. RANDOM NUMBER GENERATION: This program doesn't reseed the GAUSS random number generator for large SAMPLE you may want to occasionally reseed the generator, so that random numbers don't cycle. */ @-------------------MAIN PROGRAM----------------------------------------------@ new; cls; print off; gosub init; gosub Lmat; gosub sampleup; gosub generate; @--generates equicorrelated k-variate t--@ gosub percent; @ calculates the percentiles @ end; @-----------------------------------------------------------------------------@ @---------------------INITIALIZATION------------------------------------------@ init: k = 10; @ ultimately one variable is control so resulting @ @ simulations will be for k-1 variables @ alpha = 0.05; ctrl = 1; @ Index of the MCC control @ sample = 100000; @ Simulation sample size @ times = 2000; @ number of of rvs generated with each "generate" @ loadm varalph; @ varalph contains the general variance matrix @ print "Calculating some preliminary results"; chola = chol(varalph); @ chola is the choleski decomposition of varalph @ @ (chola)'chola = varalph @ up = sample*alpha; @ number of obs in upper alpha% of distribution @ uptail = .001*ones(up,1); @ upper alpha% of distrib. initial values = .001 @ maxabn = 0; @ holds realizations of max of |k-variate normal| @ return; @-----------------------------------------------------------------------------@ @-------------------L MATRIX GENERATION---------------------------------------@ @ this subroutine creates a matrix, L (kxk-1), So that a k vect of elements @ @ Z can be transformed into k-1 MCC variables where ctrl is the control @ @ element. L'Z = k-1 MCC variables. @ @ That is, if ctrl = k, and the parameters of interest are g'=[g(1)...g(k)] @ @ then (L'g)' = [g(k)-g(1)...g(k)-g(k-1)]. @ Lmat: count = 1; toprows = ctrl - count; botrows = k-ctrl; if toprows < 1; top = ones(1,k-1); bot = -eye(k-1); L = top|bot; elseif botrows < 1; top = -eye(k-1); bot = ones(1,k-1); L = top|bot; else; top = -eye(toprows)~zeros(toprows,botrows); mid = ones(1,toprows)~ones(1,botrows); bot = zeros(botrows,toprows)~-eye(botrows); L = top|mid|bot; endif; diagLVL = diag(L'varalph*L); clear varalph; cls; return; @----------------------------------------------------------------------------@ @-------------------------SAMPLE UP THROUGH SAMPLE SIZE----------------------@ SAMPLEUP: iter = 0; do while iter < sample; gosub generate; iter = iter + times; uptail = uptail|maxabn; @ concatonate TIMES new rvs to upper alpha% of @ @ distribution. @ rnk = rankindx(uptail,1); @ find the rank of obs in the upper alpha% @ dum = dummydn(rnk,times,1); @ a vector with 0 in the jth position if the @ @ observation in the jth position in uptail @ @ is one of the TIMES smallest observations @ @ in uptail, 1 otherwise @ clear rnk; uptail = uptail.*dum; @ replaces the TIMES smallest obs with 0 @ clear dum; uptail = delif(uptail,uptail .==0); @ removes the N smallest obs @ iter; @ prints the iteration to screen @ endo; return; @------------------------------------------------------------------------------@ @--------------------GENERATE "TIMES" REALIZATIONS-----------------------------@ generate: w = rndn(k,times); @ generates k iid standard normal variables @ @ that is w is distributed N(0, I) @ z = chola'w; @ z is distributed N(0, varalpha) @ n = L'z./sqrt(diagLVL); @ n is distibuted N(0, L'varalpaL/diag(LVL) @ maxabn = maxc(abs(n)); @ for a one-sided critical value remove "abs(.)" @ return; @-----------------------------------------------------------------------------@ @-------------------PULL CRITICAL VALUE---------------------------------------@ percent: critval = minc(uptail); print "CRITICAL VALUE max |N| IS: " critval ""; return; @-----------------------------------------------------------------------------@