function [pvalue,frac_not_pd] = pvalue_TestStat_con_GEV(TestStat,k,T);
    % Compute multiple draws for the test statistic under the null of stability and xsi = 1
    % Compute pvalue of the test statistic given by TestStat
    
    % Set up random number generator
    rng(1234);
    sigma = 1.0;
    xi = 1.0;
    mu = sigma/xi;
    N = 1000;
    TestStat_mat = zeros(N,1);

    % Set options for optimizers
    options_fminunc = optimset('Display','off');
    options_fminsearch = optimset('Display','off');

    % Generate Exponential Random Variables. These are held fixed across the value of xi
    E_mat = exprnd(1,k,T,N);
    CumSum_E_mat = cumsum(E_mat,1);  % Cumulation sum in k dimension  (This is irrelevant for k=1, but included for generality for k>1 applications)
    % Generate the GEV draws ... store in a matrix that is k x T x N               % Value of mu under the null
    z_mat = (CumSum_E_mat.^(-xi)-1)/xi;      % GEV draws with sigma = 1 and mu = 0
    y_mat = sigma*z_mat+mu;                  % GEV draws with sigma and mu
    x_start_unc = [xi sigma mu];             % Starting Values for optimizers -- use true values (unconstrained MLE)
    x_start_con = [sigma];                % Starting Values for optimizers -- use true values (constrained MLE)

    sum_not_pd = 0;
    for i_N = 1:N
        y = squeeze(y_mat(:,:,i_N));  % k x T matrix of extremes
        Y = y'; 
        
        % Compute the statistics 
        % Find MLE from GEV distribution: Order of parameters is xsi, sigma, mu 
        f_unc = @(x) minus_llf_gev(x,Y);
        [x_1,fval_1]=fminunc(f_unc,x_start_unc,options_fminunc); 
        [x_2,fval_2] = fminsearch(f_unc,x_start_unc,options_fminsearch); 
        if fval_1 < fval_2
            x_mle = x_1;
        else
            x_mle = x_2;
        end
        xsi_mle = x_mle(1); 
        sigma_mle = x_mle(2);
        mu_mle = x_mle(3);

        % Compute Hessian using unconstrained MLE
        [~,~,info_h_unc] = score_info_GEV_xi_sigma_mu(Y,xsi_mle,sigma_mle,mu_mle);

        tmp = info_h_unc;
        [info_h_unc,i_pd_flag] = make_pd(info_h_unc);  % Make sure it is positive definite
        sum_not_pd = sum_not_pd + i_pd_flag;

        % Compute MLE under the constraint xsi = 1 and mu = sigma/xi
        f_con = @(x) minus_llf_gev_con_xsi_mu(x,Y);
        [x_1,fval_1] = fminunc(f_con,x_start_con,options_fminunc);
        [x_2,fval_2] = fminsearch(f_con,x_start_con,options_fminsearch); 
        if fval_1 < fval_2
            x_mle_con = x_1;
         else
            x_mle_con = x_1;
        end
        xsi_mle_con = 1.0;
        sigma_mle_con = x_mle_con;
        mu_mle_con = sigma_mle_con/xsi_mle_con;
        
        % Compute the score under the constraint
        score_mat_con = score_GEV_xi_sigma_mu(Y,xsi_mle_con, sigma_mle_con, mu_mle_con);
          
        % Nyblom Statistic
        L_Stat_con = compute_nyblom(score_mat_con,info_h_unc);
        % LM Statistic that xsi = 1 and mu = sigma/xi
        V = info_h_unc;
        Vinv = inv(V);
        SumScore = sum(score_mat_con);
        LM_stat_con = SumScore*Vinv*SumScore'/T;
        % Combined Statistic
        TestStat_mat(i_N,1) = L_Stat_con+LM_stat_con;
    end

    % Compute p-value
    pvalue = mean(TestStat_mat > TestStat);
    frac_not_pd = sum_not_pd/N;
end

