% calculate critical values for test of regime switching
% CV based on Cho and White Asymptotic Theory
% Linear Regression Models with Gaussian Errors
% user is prompted to enter limits for the grid of values
% number of terms in the expansion is set at 150, grid mesh is set at .01
% Written by Doug Steigerwald, Spring 2011

clc;                                        % Clear Console
clear;                                      % Clear Existing Variables

d = 150;                                    % largest term in the expansion
j = 3:d;                                    % index of expansion terms

nr = 100000;                                % number of replications

% user input for the upper and lower limits for the grid of values
ul = input('Enter the upper limit for the grid of values: ');
ll = input('Enter the lower limit for the grid of values: ');

nh = floor((ul-ll)*100)+1;                  % number of terms in the grid
eta = linspace(ll,ul,nh);                   % grid of values at which to evaluate the process

stddev = rdivide(1,sqrt(gamma(j+1)))';      % vector of 1/root(j!)

emat = normrnd(0,repmat(stddev,1,nr),d-2,nr); % matrix of simulated N(0,stddev^2) random variates
                                            % # columns = nr
                                            % # rows = number of terms in
                                            % expansion (d-2)                                            
                                            
etamat = repmat(eta,d-2,1).^repmat(j',1,nh);     % create a matrix: 1 column for each value in the grid, 1 row for each term in the expansion
                                            % first column (eta(1)^3,
                                            % eta(1)^4,...,eta(1)^d) where
                                            % eta(1) is the first element
                                            % of the grid                                            

s1 = etamat'*emat;                          % matrix of sums, nh rows, nr columns (1,2): sum of eta(1)^j*e(j,rep 2) and (2,1) sum of eta(2)^j*e(j,rep 1)
  
sc = sqrt(exp(eta.^2)-1-eta.^2-eta.^4/2);   % scaling factor

cs = s1./repmat(sc',1,nr);                  % rescale matrix of sums by scaling factor
cs1 = vertcat(cs,zeros(1,nr));              % add row of zeros at end to insert the 0 for calculation of min(0,g(eta))

m2 = min(cs1).^2;                           % take minimum of each column, which is min(0,g(eta)) for each replication and construct [min(0,g(eta))]^2

m1 = max(vertcat(emat(2,1:nr)*sqrt(gamma(5)),zeros(1,nr))).^2; % construct (max(0,epslion4)^2

m = sort(max(vertcat(m1,m2)));               % max( [max(0,epsilon4)]^2, [min(0,g(eta))]^2 ), sorted in ascending order across replications
                                            
kcv = floor(.95*nr);                        % Calculate the value that corresponds to the .95 quantile of the replication vector, without bias correction

fprintf('The .95 critical value based on %g replications is \n', nr);
cv = m(kcv)                                % Construct the critical value as the .95 quantile of the sorted maxima
