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

Загрузка набора данных о потреблении электроэнергии в Канаде от Всемирного банка. Ответ - потребление электроэнергии Канадой в кВт· ч (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;

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

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

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

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

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

Коррекция ковариационной матрицы коэффициента 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, обозначенный $${}_{\mathrm{}}$$в таблице как x1, меньше обычной стандартной ошибки OLS. Это говорит о том, что в этом наборе данных коррекция остаточной гетероскедастичности и автокорреляции повышает точность измерения линейного влияния реального ВВП на потребление энергии.

Вычислите доверительные диапазоны рабочей системы (Working-Hotelling).

Вычислите 95% доверительный диапазон Working-Hoteling для каждой оценки ковариации, используя 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.

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

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 игнорирует эту картину волатильности.

Ссылки: