www.gusucode.com > econ 案例源码程序 matlab代码 > econ/PowerOfTheChowTestExample.m

    %% Power of the Chow Test
% This example shows how to estimate the power of a Chow test using a Monte
% Carlo simulation.
%% Introduction
% Statistical power is the probability of rejecting the
% null hypothesis given that it is actually false.  To estimate the power
% of a test:
% 
% # Simulate many data sets from a model that typifies the alternative
% hypothesis.
% # Test each data set.
% # Estimate the power, which is the proportion of times the test rejects
% the null hypothesis.
%
%%
% The following can compromise the power of the Chow test:
%
% * Linear model assumption departures
% * Relatively large innovation variance 
% * Using the forecast test when the sample size of the complementary
% subsample is greater than the number of coefficients in the test
% <docid:econ_ug.buxuqps-1>.
%
% Departures from model assumptions allow for an examination of the
% factors that most affect the power of the Chow test.
%%
% Consider the model
%
% $$\texttt{y} = \left[ {\begin{array}{*{20}{c}}
% {\texttt{X1}}&0\\
% 0&{\texttt{X2}}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {\texttt{beta1}}\\
% {\texttt{beta2}}
% \end{array}} \right] + \texttt{innov}$$
%
% * |innov| is a vector of random Gaussian variates with mean zero and
% standard deviation |sigma|.
% * |X1| and |X2| are the sets of predictor data for initial and
% complementary subsamples, respectively.
% * |beta1| and |beta2| are the regression coefficient vectors for the
% initial and complementary subsamples, respectively.
%
%% Simulate Predictor Data
%%
% Specify four predictors, 50 observations, and a break point at period 44
% for the simulated linear model.

% Copyright 2015 The MathWorks, Inc.

numPreds = 4;
numObs = 50;
bp = 44;
rng(1); % For reproducibility
%%
% Form the predictor data by specifying means for the predictors, and then
% adding random, standard Gaussian noise to each of the means.
mu = [0 1 2 3];
X = repmat(mu,numObs,1) + randn(numObs,numPreds);
%%
% To indicate an intercept, add a column of ones to the predictor data.
X = [ones(numObs,1) X];
X1 = X(1:bp,:);        % Initial subsample predictors
X2 = X(bp+1:end,:);    % Complementary subsample predictors
%%
% Specify the true values of the regression coefficients.
beta1 = [1 2 3 4 5]'; % Initial subsample coefficients
%% Estimate Power for Small and Large Jump
% Compare the power between the break point and forecast tests for jumps of
% different sizes small in the intercept and second regression coefficient.
% In this example, a small jump is a 10% increase in the current value, and
% a large jump is a 15% increase. Complementary subsample coefficients
beta2Small = beta1 + [beta1(1)*0.1 0 beta1(3)*0.1 0 0 ]'; 
beta2Large = beta1 + [beta1(1)*0.15 0 beta1(3)*0.15 0 0 ]'; 
%%
% Simulate 1000 response paths of the linear model for each of the small
% and large coefficient jumps. Specify that |sigma| is 0.2.  Choose to test
% the intercept and the second regression coefficient.
M = 1000;

sigma = 0.2;
Coeffs = [true false true false false];
h1BP = nan(M,2); % Preallocation
h1F = nan(M,2);
for j = 1:M
    innovSmall = sigma*randn(numObs,1);
    innovLarge = sigma*randn(numObs,1);
    ySmall = [X1 zeros(bp,size(X2,2)); ...
        zeros(numObs - bp,size(X1,2)) X2]*[beta1; beta2Small] + innovSmall;
    yLarge = [X1 zeros(bp,size(X2,2)); ...
        zeros(numObs - bp,size(X1,2)) X2]*[beta1; beta2Large] + innovLarge;
    h1BP(j,1) = chowtest(X,ySmall,bp,'Intercept',false,'Coeffs',Coeffs,...
        'Display','off')';
    h1BP(j,2) = chowtest(X,yLarge,bp,'Intercept',false,'Coeffs',Coeffs,...
        'Display','off')';
    h1F(j,1) = chowtest(X,ySmall,bp,'Intercept',false,'Coeffs',Coeffs,...
        'Test','forecast','Display','off')';
    h1F(j,2) = chowtest(X,yLarge,bp,'Intercept',false,'Coeffs',Coeffs,...
        'Test','forecast','Display','off')';
end
%%
% Estimate the power by computing the proportion of times |chowtest|
% correctly rejected the null hypothesis of coefficient stability.
power1BP = mean(h1BP);
power1F = mean(h1F);
table(power1BP',power1F','RowNames',{'Small_Jump','Large_Jump'},...
    'VariableNames',{'Breakpoint','Forecast'})
%%
% In this scenario, the Chow test can detect a change in the coefficient
% with more power when the jump is larger.  The break point test has
% greater power to detect the jump than the forecast test.
%% Estimate Power for Large Innovations Variance
% Simulate 1000 response paths of the linear model for a large coefficient
% jump. Specify that |sigma| is 0.4.  Choose to test the intercept and the
% second regression coefficient.
sigma = 0.4;
h2BP = nan(M,1); 
h2F = nan(M,1);
for j = 1:M
    innov = sigma*randn(numObs,1);
    y = [X1 zeros(bp,size(X2,2)); ...
        zeros(numObs - bp,size(X1,2)) X2]*[beta1; beta2Large] + innov;
    h2BP(j) = chowtest(X,y,bp,'Intercept',false,'Coeffs',Coeffs,...
        'Display','off')';
    h2F(j) = chowtest(X,y,bp,'Intercept',false,'Coeffs',Coeffs,...
        'Test','forecast','Display','off')';
end

power2BP = mean(h2BP);
power2F = mean(h2F);
table([power1BP(2); power2BP],[power1F(2); power2F],...
    'RowNames',{'Small_sigma','Large_Sigma'},...
    'VariableNames',{'Breakpoint','Forecast'})
%%
% For larger innovation variance, both Chow tests have difficulty detecting
% the large structural breaks in the intercept and second regression
% coefficient.