Обзор тематического исследования

Этот раздел содержит пример рабочего процесса, описанного в Рабочем процессе Модели VAR. Пример использует три временных рядов: GDP, денежная масса M1 и 3-месячный уровень Казначейского векселя. Пример показывает:

Загрузите и преобразуйте данные

Файл Data_USEconModel поставки с программным обеспечением Econometrics Toolbox™. Файл содержит временные ряды от базы данных Federal Reserve Bank of St. Louis Economics Data (FRED) в табличном массиве. Этот пример использует три из временных рядов:

Загрузите набор данных. Создайте переменную для действительного GDP.

load Data_USEconModel
DataTable.RGDP = DataTable.GDP./DataTable.GDPDEF*100;

Отобразите данные на графике, чтобы искать тренды.

figure
subplot(3,1,1)
plot(DataTable.Time,DataTable.RGDP,'r');
title('Real GDP')
grid on
subplot(3,1,2);
plot(DataTable.Time,DataTable.M1SL,'b');
title('M1')
grid on
subplot(3,1,3);
plot(DataTable.Time,DataTable.TB3MS,'k')
title('3-mo T-bill')
grid on

Действительный GDP и данные M1, кажется, растут экспоненциально, в то время как Казначейский вексель возвращается, не показывают экспоненциального роста. Чтобы противостоять трендам в действительном GDP и M1, возьмите различие логарифмов данных. Кроме того, стабилизируйте ряд Казначейского векселя путем взятия первого различия. Синхронизируйте ряд даты так, чтобы данные имели одинаковое число строк для каждого столбца.

rgdpg = price2ret(DataTable.RGDP);
m1slg = price2ret(DataTable.M1SL);
dtb3ms = diff(DataTable.TB3MS);
Data = array2timetable([rgdpg m1slg dtb3ms],...
    'RowTimes',DataTable.Time(2:end),'VariableNames',{'RGDP' 'M1SL' 'TB3MS'});

figure
subplot(3,1,1)
plot(Data.Time,Data.RGDP,'r');
title('Real GDP')
grid on
subplot(3,1,2);
plot(Data.Time,Data.M1SL,'b');
title('M1')
grid on
subplot(3,1,3);
plot(Data.Time,Data.TB3MS,'k'),
title('3-mo T-bill')
grid on

Шкала первых двух столбцов приблизительно в 100 раз меньше, чем третье. Умножьте первые два столбца на 100 так, чтобы временные ряды были всеми примерно по той же шкале. Это масштабирование дает возможность строить весь ряд на том же графике. Что еще более важно, этот тип масштабирования делает оптимизацию более численно устойчивой (например, максимизируя логарифмическую правдоподобность).

Data{:,1:2} = 100*Data{:,1:2};
figure
plot(Data.Time,Data.RGDP,'r');
hold on
plot(Data.Time,Data.M1SL,'b'); 
datetick('x') 
grid on
plot(Data.Time,Data.TB3MS,'k');
legend('Real GDP','M1','3-mo T-bill'); 
hold off

Выберите и подбирайте модели

Можно выбрать много различных моделей для данных. Этот пример использует четыре модели.

Удалите все отсутствующие значения с начала ряда.

idx = all(~ismissing(Data),2);
Data = Data(idx,:);

Создайте эти четыре модели.

numseries = 3;
dnan = diag(nan(numseries,1));
seriesnames = {'Real GDP','M1','3-mo T-bill'};
VAR2diag = varm('AR',{dnan dnan},'SeriesNames',seriesnames);
VAR2full = varm(numseries,2);
VAR2full.SeriesNames = seriesnames;
VAR4diag = varm('AR',{dnan dnan dnan dnan},'SeriesNames',seriesnames);
VAR4full = varm(numseries,4);
VAR4full.SeriesNames = seriesnames;

Матричный dnan диагональная матрица с NaN значения по его основной диагонали. В общем случае отсутствующие значения задают присутствие параметра в модели и указывают, что параметр должен быть подходящим к данным. MATLAB® содержит от диагональных элементов, 0, зафиксированный во время оценки. В отличие от этого технические требования для VAR2full и VAR4full составьте матрицы NaN значения. Поэтому estimate соответствует полным матрицам для авторегрессивных матриц.

Чтобы оценить качество моделей, создайте векторы индекса, которые делят данные об ответе на три периода: предварительная выборка, оценка и прогноз. Подбирайте модели к данным об оценке, использование преддемонстрационного периода, чтобы обеспечить изолировало данные. Сравните предсказания подобранных моделей к данным о прогнозе. Период оценки находится в выборке, и период прогноза вне выборки (также известен backtesting).

Для двух моделей VAR (4) преддемонстрационный период является первыми четырьмя строками Data. Используйте тот же преддемонстрационный период для моделей VAR (2) так, чтобы все модели были подходящими к тем же данным. Это необходимо для подходящих сравнений модели. Для обеих моделей период прогноза составляет итоговые 10% строк Data. Период оценки для моделей идет от строки 5 до 90%-й строки. Задайте эти периоды данных.

idxPre = 1:4;
T = ceil(.9*size(Data,1));
idxEst = 5:T;
idxF = (T+1):size(Data,1);
fh = numel(idxF);

Теперь, когда модели и временные ряды существуют, можно легко подбирать модели к данным.

[EstMdl1,EstSE1,logL1,E1] = estimate(VAR2diag,Data{idxEst,:},...
    'Y0',Data{idxPre,:});
[EstMdl2,EstSE2,logL2,E2] = estimate(VAR2full,Data{idxEst,:},...
    'Y0',Data{idxPre,:});
[EstMdl3,EstSE3,logL3,E3] = estimate(VAR4diag,Data{idxEst,:},...
    'Y0',Data{idxPre,:});
[EstMdl4,EstSE4,logL4,E4] = estimate(VAR4full,Data{idxEst,:},...
    'Y0',Data{idxPre,:});

Проверяйте соответствие модели

Можно проверять, являются ли предполагаемые модели устойчивыми и обратимыми путем отображения Description свойство каждого объекта. (В этих моделях нет никаких терминов MA, таким образом, модели являются обязательно обратимыми.) Описания показывают, что все предполагаемые модели устойчивы.

EstMdl1.Description

ans = 
"AR-Stationary 3-Dimensional VAR(2) Model"

EstMdl2.Description

ans = 
"AR-Stationary 3-Dimensional VAR(2) Model"

EstMdl3.Description

ans = 
"AR-Stationary 3-Dimensional VAR(4) Model"

EstMdl4.Description

ans = 
"AR-Stationary 3-Dimensional VAR(4) Model"

AR-Stationary появляется в выходе, указывающем, что авторегрессивные процессы устойчивы.

Можно сравнить ограниченные (диагональные) модели AR с их неограниченными (полными) дубликатами с помощью lratiotest. Тестовые отклонения или сбои, чтобы отклонить гипотезу, что ограниченные модели соответствуют с 5%-м допуском по умолчанию. Это - тест в выборке.

Примените тесты отношения правдоподобия. Необходимо извлечь количество предполагаемых параметров от итоговой структуры, возвращенной summarize. Затем передайте различия в количестве предполагаемых параметров и логарифмической правдоподобности к lratiotest выполнять тесты.

results1 = summarize(EstMdl1);
np1 = results1.NumEstimatedParameters;
results2 = summarize(EstMdl2);
np2 = results2.NumEstimatedParameters;
results3 = summarize(EstMdl3);
np3 = results3.NumEstimatedParameters;
results4 = summarize(EstMdl4);
np4 = results4.NumEstimatedParameters;

reject1 = lratiotest(logL2,logL1,np2 - np1)

reject1 = logical
   1

reject3 = lratiotest(logL4,logL3,np4 - np3)

reject3 = logical
   1

reject4 = lratiotest(logL4,logL2,np4 - np2)

reject4 = logical
   0

1 результаты показывают, что тест отношения правдоподобия отклонил обоих ограниченные модели в пользу соответствующих неограниченных моделей. Поэтому на основе этого теста, неограниченные модели VAR (2) и VAR (4) предпочтительны. Однако тест не отклоняет неограниченную модель VAR (2) в пользу неограниченной модели VAR (4). (Этот тест рассматривает модель VAR (2) как модель VAR (4) с ограничениями, что AR матриц авторегрессии (3) и AR (4) 0.) Поэтому кажется, что неограниченная модель VAR (2) является лучшей моделью.

Чтобы найти лучшую модель в наборе, минимизируйте Критерий информации о Akaike (AIC). Используйте в выборочных данных, чтобы вычислить AIC. Вычислите критерий этих четырех моделей.

AIC = aicbic([logL1 logL2 logL3 logL4],[np1 np2 np3 np4])

AIC = 1×4
103 ×

    1.4794    1.4396    1.4785    1.4537

Лучшая модель согласно этому критерию является неограниченной моделью VAR (2). Заметьте также, что неограниченная модель VAR (4) имеет более низкую информацию о Akaike, чем любая из ограниченных моделей. На основе этого критерия неограниченная модель VAR (2) является лучшей с неограниченной моделью VAR (4), появляющейся затем в настройку.

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

[FY1,FYCov1] = forecast(EstMdl1,fh,Data{idxEst,:});
[FY2,FYCov2] = forecast(EstMdl2,fh,Data{idxEst,:});
[FY3,FYCov3] = forecast(EstMdl3,fh,Data{idxEst,:});
[FY4,FYCov4] = forecast(EstMdl4,fh,Data{idxEst,:});

Оцените аппроксимированные 95%-е интервалы прогноза для модели оптимальной подгонки.

extractMSE = @(x)diag(x)';
MSE = cellfun(extractMSE,FYCov2,'UniformOutput',false);
SE = sqrt(cell2mat(MSE));

YFI = zeros(fh,EstMdl2.NumSeries,2);
YFI(:,:,1) = FY2 - 2*SE;
YFI(:,:,2) = FY2 + 2*SE;

Этот график показывает предсказания модели оптимальной подгонки в теневой области направо.

figure;
for j = 1:EstMdl2.NumSeries
    subplot(3,1,j);
    h1 = plot(Data.Time((end-49):end),Data{(end-49):end,j});
    hold on;
    h2 = plot(Data.Time(idxF),FY2(:,j));
    h3 = plot(Data.Time(idxF),YFI(:,j,1),'k--');
    plot(Data.Time(idxF),YFI(:,j,2),'k--');
    title(EstMdl2.SeriesNames{j});
    h = gca;
    fill([Data.Time(idxF(1)) h.XLim([2 2]) Data.Time(idxF(1))],...
        h.YLim([1 1 2 2]),'k','FaceAlpha',0.1,'EdgeColor','none');
    legend([h1 h2 h3],'True','Forecast','95% Forecast interval',...
        'Location','northwest')
    hold off;
end

Это теперь прямо, чтобы вычислить ошибку суммы квадратов между предсказаниями и данными.

Error1 = Data{idxF,:} - FY1;
Error2 = Data{idxF,:} - FY2;
Error3 = Data{idxF,:} - FY3;
Error4 = Data{idxF,:} - FY4;

SSerror1 = Error1(:)' * Error1(:);
SSerror2 = Error2(:)' * Error2(:);
SSerror3 = Error3(:)' * Error3(:);
SSerror4 = Error4(:)' * Error4(:);
figure
bar([SSerror1 SSerror2 SSerror3 SSerror4],.5)
ylabel('Sum of squared errors')
set(gca,'XTickLabel',...
    {'AR2 diag' 'AR2 full' 'AR4 diag' 'AR4 full'})
title('Sum of Squared Forecast Errors')

Прогнозирующая эффективность этих четырех моделей подобна.

Полная модель AR (2), кажется, является лучшей и большая часть экономной подгонки. Его параметры модели следующие.

summarize(EstMdl2)

 
   AR-Stationary 3-Dimensional VAR(2) Model
 
    Effective Sample Size: 176
    Number of Estimated Parameters: 21
    LogLikelihood: -698.801
    AIC: 1439.6
    BIC: 1506.18
 
                     Value       StandardError    TStatistic      PValue  
                   __________    _____________    __________    __________

    Constant(1)       0.34832       0.11527          3.0217      0.0025132
    Constant(2)       0.55838        0.1488          3.7526     0.00017502
    Constant(3)      -0.45434       0.15245         -2.9803      0.0028793
    AR{1}(1,1)        0.26252       0.07397          3.5491     0.00038661
    AR{1}(2,1)      -0.029371      0.095485         -0.3076        0.75839
    AR{1}(3,1)        0.22324      0.097824          2.2821       0.022484
    AR{1}(1,2)      -0.074627      0.054476         -1.3699        0.17071
    AR{1}(2,2)         0.2531      0.070321          3.5992     0.00031915
    AR{1}(3,2)      -0.017245      0.072044        -0.23936        0.81082
    AR{1}(1,3)       0.032692      0.056182         0.58189        0.56064
    AR{1}(2,3)       -0.35827      0.072523           -4.94     7.8112e-07
    AR{1}(3,3)       -0.29179        0.0743         -3.9272     8.5943e-05
    AR{2}(1,1)        0.21378      0.071283          2.9991      0.0027081
    AR{2}(2,1)      -0.078493      0.092016        -0.85304        0.39364
    AR{2}(3,1)        0.24919      0.094271          2.6433      0.0082093
    AR{2}(1,2)        0.13137      0.051691          2.5415       0.011038
    AR{2}(2,2)        0.38189      0.066726          5.7233      1.045e-08
    AR{2}(3,2)       0.049403      0.068361         0.72269        0.46987
    AR{2}(1,3)       -0.22794      0.059203           -3.85     0.00011809
    AR{2}(2,3)     -0.0052932      0.076423       -0.069262        0.94478
    AR{2}(3,3)       -0.37109      0.078296         -4.7397     2.1408e-06

 
   Innovations Covariance Matrix:
    0.5931    0.0611    0.1705
    0.0611    0.9882   -0.1217
    0.1705   -0.1217    1.0372

 
   Innovations Correlation Matrix:
    1.0000    0.0798    0.2174
    0.0798    1.0000   -0.1202
    0.2174   -0.1202    1.0000

Предскажите наблюдения

Можно сделать предсказания или прогнозы с помощью подобранной модели (EstMdl2) любой:

В обоих случаях преобразуйте прогнозы, таким образом, они непосредственно сопоставимы с исходными временными рядами.

Сгенерируйте 10 предсказаний от подобранной модели, начинающейся в последние времена с помощью forecast.

[YPred,YCov] = forecast(EstMdl2,10,Data{idxF,:});

Преобразуйте предсказания путем отмены масштабирования, и дифференцирование применилось к исходным данным. Убедитесь, что вставили последнее наблюдение в начале временных рядов перед использованием cumsum отменить дифференцирование. И, поскольку дифференцирование произошло после взятия логарифмов, вставьте логарифм перед использованием cumsum.

YFirst = DataTable(idx,{'RGDP' 'M1SL' 'TB3MS'}); 
EndPt = YFirst{end,:};
EndPt(:,1:2) = log(EndPt(:,1:2));
YPred(:,1:2) = YPred(:,1:2)/100;     % Rescale percentage
YPred = [EndPt; YPred];              % Prepare for cumsum
YPred(:,1:3) = cumsum(YPred(:,1:3));
YPred(:,1:2) = exp(YPred(:,1:2));
fdates = dateshift(YFirst.Time(end),'end','quarter',0:10); % Insert forecast horizon

figure
for j = 1:EstMdl2.NumSeries
    subplot(3,1,j)
    plot(fdates,YPred(:,j),'--b')
    hold on
    plot(YFirst.Time,YFirst{:,j},'k')
    grid on 
    title(EstMdl2.SeriesNames{j})
    h = gca;
    fill([fdates(1) h.XLim([2 2]) fdates(1)],h.YLim([1 1 2 2]),'k',...
        'FaceAlpha',0.1,'EdgeColor','none');
    hold off
end

Графики показывают экстраполяции синими пунктирными линиями в светло-сером горизонте прогноза и исходным рядом данных в чистом черном.

Посмотрите на последние несколько лет в этом графике получить смысл того, как предсказания относятся к последним точкам данных.

YLast = YFirst(170:end,:);

figure
for j = 1:EstMdl2.NumSeries
    subplot(3,1,j)
    plot(fdates,YPred(:,j),'b--')
    hold on
    plot(YLast.Time,YLast{:,j},'k')
    grid on
    title(EstMdl2.SeriesNames{j})
    h = gca;
    fill([fdates(1) h.XLim([2 2]) fdates(1)],h.YLim([1 1 2 2]),'k',...
        'FaceAlpha',0.1,'EdgeColor','none');
    hold off
end

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

В качестве альтернативы можно сгенерировать 10 предсказаний от подобранной модели, начинающейся в последние времена с помощью simulate. Этот метод симулирует 2000 раз временных рядов, и затем генерирует средние значения и стандартные отклонения в течение каждого периода. Средние значения отклонения в течение каждого периода являются предсказаниями в течение того периода.

Симулируйте временные ряды от подобранной модели, начинающейся в последние времена.

rng(1); % For reproducibility
YSim = simulate(EstMdl2,10,'Y0',Data{idxF,:},'NumPaths',2000);

Преобразуйте предсказания путем отмены масштабирования, и дифференцирование применилось к исходным данным. Убедитесь, что вставили последнее наблюдение в начале временных рядов перед использованием cumsum отменить дифференцирование. И, поскольку дифференцирование произошло после взятия логарифмов, вставьте логарифм перед использованием cumsum.

EndPt = YFirst{end,:}; 
EndPt(1:2) = log(EndPt(1:2)); 
YSim(:,1:2,:) = YSim(:,1:2,:)/100; 
YSim = [repmat(EndPt,[1,1,2000]);YSim]; 
YSim(:,1:3,:) = cumsum(YSim(:,1:3,:)); 
YSim(:,1:2,:) = exp(YSim(:,1:2,:));

Вычислите среднее и стандартное отклонение каждого ряда и постройте результаты. График имеет среднее значение черного цвета цвета, с + стандартное отклонение/-1 красного цвета.

YMean = mean(YSim,3);
YSTD = std(YSim,0,3);

figure
for j = 1:EstMdl2.NumSeries
    subplot(3,1,j)
    plot(fdates,YMean(:,j),'k')
    grid on
    hold on
    plot(YLast.Time(end-10:end),YLast{end-10:end,j},'k')
    plot(fdates,YMean(:,j) + YSTD(:,j),'--r')
    plot(fdates,YMean(:,j) - YSTD(:,j),'--r')
    title(EstMdl2.SeriesNames{j})
    h = gca;
    fill([fdates(1) h.XLim([2 2]) fdates(1)],h.YLim([1 1 2 2]),'k',...
        'FaceAlpha',0.1,'EdgeColor','none');
    hold off
end

Графики показывают увеличивающийся рост в GDP, умеренном к небольшому росту в M1 и неопределенности по поводу направления уровней Казначейского векселя.