%% 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 [24].
% 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.

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.
% Copyright 2015 The MathWorks, Inc.