Использование bootstrap в Matlab для проведения вероятного теста отношения и вычисления значения p и мощности теста - PullRequest
Новая
       29

Использование bootstrap в Matlab для проведения вероятного теста отношения и вычисления значения p и мощности теста

0 голосов
/ 23 марта 2020

Я новичок в Matlab. Сначала я имитирую данные 1000 раз и запускаю регрессию по X и Y (как OLS, так и ограниченным OLS). Для каждого смоделированного набора данных я генерирую p-значение для отношения правдоподобия c и вычисляю мощность теста. Я получил правильный ответ для этой части. В части bootstrap сначала выполняется процедура получения остатка из последних смоделированных данных. Затем выполните повторную выборку остатка 1000 раз и добавьте новые остатки в X * betahat (который получен из исходной регрессии OLS), чтобы получить новую переменную Yb. Затем регрессируйте Yb и переменную X (как OLS, так и ограниченную OLS), чтобы получить статистику отношения правдоподобия c для каждого bootstrap. Однако в этой части я получаю очень большое значение p, и тестовая мощность теста LR близка к нулю.

Первая часть - это симуляция, в которой я получаю правильный ответ. 

randn('state',7)
N = 1000;
sig = 3;
T = 40;
cc = T/8;
beta = [0 -5 3 5]';
dum1 = [ones(T/2,1);zeros(T/2,1)]; dum2 = 1 - dum1;
time = kron((1:4)',ones(cc,1));
c3 = kron([1,0]',time); c4 = kron([0,1]',time);
X = [dum1 dum2 c3 c4];
H = [1 -1 0 0; 0 0 1 -1];
A = inv(X'*X);
Ftestresults = zeros(N,1);
LRtestresults = zeros(N,1);
Z = [dum1+dum2,c3+c4];
for i = 1:N
    Y = X*beta + sig*randn(T,1); 
    betahat = A*X'*Y;
    Sbeta = (Y-X*betahat)'*(Y-X*betahat);
    gammahat = OLSrestrict(Y,X,H);
    yhat = X*gammahat;
    res = Y - yhat;
    Sgamma = sum(res.^2);
    sigmahat = Sbeta/(T-4);
    Fvalue = (Sgamma-Sbeta)/(2*sigmahat);
    pvalue = 1-fcdf(Fvalue,2,T-4);
    LR = -2*log((Sgamma/Sbeta)^(-T/2));
    pvalue2 = 1-chi2cdf(LR,2);
    if pvalue < 0.05
        Ftestresults(i) = 1;
    end
    if pvalue2 < 0.05
        LRtestresults(i) = 1;
    end
end
Fpower = sum(Ftestresults)/N
LRpower = sum(LRtestresults)/N
function gamma = OLSrestrict(y,X,H,h)
    [J,K] = size(H);
    if nargin <4,
        h = zeros(J,1);
    end
    b = regress(y,X);
    A = inv(X'*X);
    gamma = b + A*H'*inv(H*A*H')*(h-H*b);
end

Вторая часть - bootstrap. 

function pvaluebootstrap = bootstrappvalue(Y,X,T,H,LR)
    betahat = inv(X'*X)*X'*Y;
    residual = Y - X*betahat;
    B = 1000;
    LRdis = zeros(B,1);
    for b = 1:B
       i = randi([1,T],1,T);
       residualb = zeros(T,1);
       for j = 1:T
           residualb(j) = residual(i(j));
       end
       Yb = X*betahat+residualb;
       betanew = inv(X'*X)*X'*Yb;
       res1 = Yb - X*betanew;
       Sbetanew = sum(res1.^2);
       gammanew = OLSrestrict(Yb,X,H); 
       res2 = Yb-X*gammanew; Sgammanew = sum(res2.^2);
       LRdis(b) = -2*log((Sgammanew/Sbetanew)^(-T/2));
    end
    count = 0;
    for k = 1:B
        if LRdis(k)>LR
            count = count+1;
        end
    end
    pvaluebootstrap = count/B;
end

Я действительно смущен и не могу понять, что происходит не так. Спасибо за помощь!

...