Загрузите данные

Загрузите канадский набор данных потребления электроэнергии от World Bank. Ответ является потреблением электроэнергии Канады в kWh (DataTable.consump), предиктор является GDP Канады в году 2 000 долларов США (DataTable.gdp), и набор данных также содержит дефлятор GDP (DataTable.gdp_deflator). Поскольку DataTable расписание, DataTable.Time демонстрационный год.

load Data_PowerConsumption

Задайте модель

Смоделируйте поведение ежегодного различия в потреблении электроэнергии относительно действительного GDP как линейная модель:

consumpDiff = DataTable.consump - lagmatrix(DataTable.consump,1); ...
    % Annual difference in consumption
T = size(consumpDiff,1);
consumpDiff = consumpDiff(2:end)/1.0e+10; 
    % Scale for numerical stability
rGDP = DataTable.gdp./(DataTable.gdp_deflator);    % Deflate GDP
rGDP = rGDP(2:end)/1.0e+10;   % Scale for numerical stability
rGDPdes = [ones(T-1,1) rGDP]; % Design matrix

Mdl = fitlm(rGDP,consumpDiff);
coeff = Mdl.Coefficients(:,1);
EstParamCov = Mdl.CoefficientCovariance;
resid = Mdl.Residuals.Raw;

Отобразите данные на графике

Постройте различие в потреблении энергии, consumpDiff по сравнению с действительным GDP, чтобы проверять на возможный heteroscedasticity.

figure
plot(rGDP,consumpDiff,'.')
title 'Annual Difference in Energy Consumption vs real GDP - Canada';
xlabel 'real GDP (year 2000 USD)';
ylabel 'Annual Difference in Energy Consumption (kWh)';

Фигура указывает, что heteroscedasticity сила присутствует в ежегодном различии в потреблении энергии. Когда действительный GDP увеличивается, ежегодное различие в потреблении энергии, кажется, меньше переменной.

Постройте остаточные значения.

Постройте остаточные значения Mdl против подходящих значений и год, чтобы оценить heteroscedasticity и автокорреляцию.

figure
subplot(2,1,1)
hold on
plot(Mdl.Fitted,resid,'.')
plot([min(Mdl.Fitted) max(Mdl.Fitted)],[0 0],'k-')
title 'Residual Plots';
xlabel 'Fitted Consumption';
ylabel 'Residuals';
axis tight
hold off
subplot(2,2,3)
autocorr(resid)
h1 = gca;
h1.FontSize = 8;
subplot(2,2,4)
parcorr(resid)
h2 = gca;
h2.FontSize = 8;

Остаточный график показывает уменьшающееся остаточное отклонение с увеличением подходящего потребления. Автокорреляционная функция показывает, что автокорреляция может присутствовать в первых нескольких изолированных остаточных значениях.

Протестируйте на heteroscedasticity и автокорреляцию.

Протестируйте на условное выражение heteroscedasticity использование теста ДУГИ Энгла. Протестируйте на автокорреляцию с помощью теста Ljung-поля Q. Протестируйте на полную корреляцию с помощью теста Дербин-Уотсона.

[~,englePValue] = archtest(resid);
englePValue

englePValue = 0.1463

[~,lbqPValue] = lbqtest(resid,'lags',1:3);...
    % Significance of first three lags
lbqPValue

lbqPValue = 1×3

    0.0905    0.1966    0.0522

[dwPValue] = dwtest(Mdl);
dwPValue

dwPValue = 0.0024

P значение теста ДУГИ Энгла предлагает значительное условное выражение heteroscedasticity на 15%-м уровне значения. P значение для теста Ljung-поля Q предлагает значительную автокорреляцию с первыми и третьими изолированными остаточными значениями на 10%-м уровне значения. P значение для теста Дербин-Уотсона предполагает, что существуют убедительные доказательства для полной остаточной автокорреляции. Результаты тестов предполагают, что стандартные линейные условия модели homoscedasticity и некоррелированых ошибок нарушены, и выводы на основе содействующей ковариационной матрицы OLS являются подозреваемым.

Один способ возобновить вывод (такой как построение полосы доверия) состоит в том, чтобы откорректировать содействующую ковариационную матрицу OLS путем оценки Newey-западной содействующей ковариации.

Оцените Newey-западную содействующую ковариацию.

Откорректируйте содействующую ковариационную матрицу OLS путем оценки Newey-западной содействующей ковариации с помощью hac. Вычислите максимальную задержку, которая будет взвешена для стандартной Newey-западной оценки, maxLag (Newey и West, 1994). Используйте hac оценить стандартную Newey-западную содействующую ковариацию.

maxLag = floor(4*(T/100)^(2/9));
[NWEstParamCov,~,NWCoeff] = hac(Mdl,'type','hac',...
    'bandwidth',maxLag + 1);

Estimator type: HAC
Estimation method: BT
Bandwidth: 4.0000
Whitening order: 0
Effective sample size: 49
Small sample correction: on

Coefficient Covariances:

       |  Const      x1   
--------------------------
 Const |  0.3720  -0.2990 
 x1    | -0.2990   0.2454 

Newey-западная стандартная погрешность для коэффициента rGDP, помеченный $$x_1$$ в таблице, меньше обычной стандартной погрешности OLS. Это предполагает, что, в этом наборе данных, корректирующем для невязки heteroscedasticity и автокорреляции, увеличивает точность в измерении линейного эффекта действительного GDP на потреблении энергии.

Вычислите работающие-Hotelling полосы доверия.

Вычислите 95%, работающих-Hotelling полоса доверия на каждую оценку ковариации с помощью nlpredci (Kutner и др., 2005).

modelfun = @(b,x)(b(1)*x(:,1)+b(2)*x(:,2));
    % Define the linear model
[beta,nlresid,~,EstParamCov] = nlinfit(rGDPdes,...
    consumpDiff,modelfun,[1,1]); % estimate the model
[fity,fitcb] = nlpredci(modelfun,rGDPdes,beta,nlresid,...
    'Covar',EstParamCov,'SimOpt','on');
    % Margin of errors
conbandnl = [fity - fitcb fity + fitcb]; 
    % Confidence bands
[fity,NWfitcb] = nlpredci(modelfun,rGDPdes,...
    beta,nlresid,'Covar',NWEstParamCov,'SimOpt','on');
    % Corrected margin of error
NWconbandnl = [fity - NWfitcb fity + NWfitcb];
    % Corrected confidence bands

Постройте работающие-Hotelling полосы доверия.

Постройте работающие-Hotelling полосы доверия на тех же осях дважды: один график, отображающий потребление электроэнергии относительно действительного GDP и другое отображение временных рядов потребления электроэнергии.

figure
hold on
l1 = plot(rGDP,consumpDiff,'k.');
l2 = plot(rGDP,fity,'b-','LineWidth',2);
l3 = plot(rGDP,conbandnl,'r-');
l4 = plot(rGDP,NWconbandnl,'g--');
title 'Data with 95% Working-Hotelling Conf. Bands';
xlabel 'real GDP (year 2000 USD)';
ylabel 'Consumption (kWh)';
axis([0.7 1.4 -2 2.5])
legend([l1 l2 l3(1) l4(1)],'Data','Fitted','95% Conf. Band',...
    'Newey-West 95% Conf. Band','Location','SouthEast')
hold off

figure
year = DataTable.Time(2:end);
hold on
l1 = plot(year,consumpDiff);
l2 = plot(year,fity,'k-','LineWidth',2);
l3 = plot(year,conbandnl,'r-');
l4 = plot(year,NWconbandnl,'g--');
title 'Consumption with 95% Working-Hotelling Conf. Bands';
xlabel 'Year';
ylabel 'Consumption (kWh)';
legend([l1 l2 l3(1) l4(1)],'Consumption','Fitted',...
    '95% Conf. Band','Newey-West 95% Conf. Band',...
    'Location','SouthWest')
hold off

Графики показывают, что Newey-западное средство оценки составляет heteroscedasticity в этом, полоса доверия является широкой в областях высокой волатильности и тонкой в областях низкой волатильности. Содействующее средство оценки ковариации OLS игнорирует этот шаблон энергозависимости.

Ссылки: