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.