Загрузка данных

Загрузите канадский набор данных о потреблении электрической степени из World Bank. Ответом является потребление электроэнергии в Канаде в кВтч (DataTable.consump), предиктором является ВВП Канады в 2000 году USD (DataTable.gdp), и набор данных также содержит дефлятор ВВП (DataTable.gdp_deflator). Потому что DataTable является расписанием, DataTable.Time является годом расчета.

load Data_PowerConsumption

Определите модель

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

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 по сравнению с реальным ВВП, чтобы проверить на возможную гетероскедастичность.

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)';

Рисунок показывает, что гетероскедастичность может присутствовать в годовом различии в потреблении энергии. Когда реальный ВВП увеличивается, годовое различие в потреблении энергии, по-видимому, менее изменчива.

Постройте график невязок.

Постройте график невязок из Mdl по подобранным значениям и году для оценки гетероскедастичности и автокорреляции.

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;

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

Тест на гетероскедастичность и автокорреляцию.

Тест на условную гетероскедастичность с использованием теста ARCH Engle. Тест на автокорреляцию с использованием Q-теста Ljung-Box. Тест на общую корреляцию с использованием теста Дурбина-Ватсона.

[~,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 теста ARCH Engle предполагает значительную условную гетероскедастичность при уровне значимости 15%. Значение p для Q-теста Ljung-Box предполагает значительную автокорреляцию с первыми и третьими отстающими невязками при уровне значимости 10%. Значение p для теста Дурбина-Ватсона предполагает, что существуют убедительные доказательства общей остаточной автокорреляции. Результаты тестов предполагают, что стандартные линейные модельные условия гомосцедастичности и некоррелированные ошибки нарушены, и выводы, основанные на ковариационной матрице коэффициента OLS, подозрительны.

Один из способов продолжить вывод (такой как построение доверия диапазона) - исправить ковариационную матрицу коэффициента OLS путем оценки ковариации коэффициента Ньюи-Запада.

Оцените ковариацию коэффициента Ньюи-Уэста.

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

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 

Стандартная ошибка Ньюи-Уэста для коэффициента rGDP, маркированный $$x_1$$ в таблице меньше обычной стандартной ошибки OLS. Это говорит о том, что в этом наборе данных коррекция остаточной гетероскедастичности и автокорреляции повышает точность измерения линейного эффекта реального ВВП на энергопотребление.

Вычислите доверительные полосы Рабочей функции.

Вычислите диапазон 95% Рабочей-Hotelling доверия для каждой оценки ковариации с помощью nlpredci (Kutner et al., 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

Постройте график доверительных полос Working-Hotelling.

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

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

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

Ссылки: