Тематическое исследование модели VAR

Этот пример показывает, как анализировать модель VAR.

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

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

Загрузка и преобразование данных для стационарности
Деля преобразованные данные в предварительную выборку, оценку и интервалы прогноза, чтобы поддержать эксперимент backtesting
Создание нескольких моделей
Подбор кривой моделям к данным
Решение, какая из моделей является лучшей
Создание прогнозов на основе лучшей модели

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

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

GDP (GDP)
Денежная масса M1 (M1SL)
3-месячный уровень Казначейского векселя (TB3MS)

Загрузите набор данных. Создайте переменную для действительного 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 так, чтобы временные ряды были всеми примерно в той же шкале. Это масштабирование дает возможность строить весь ряд на том же графике. Что еще более важно, этот тип масштабирования делает оптимизацию более численно стабильной (например, максимизируя loglikelihoods).

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

Выберите и соответствуйте моделям

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

VAR (2) с авторегрессивной диагональю
VAR (2) с авторегрессивным полным
VAR (4) с авторегрессивной диагональю
VAR (4) с авторегрессивным полным

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

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,:});

Объекты модели EstMdl являются подобранными моделями.
Структуры EstSE содержат стандартные погрешности подобранных моделей.
Значения logL являются loglikelihoods подобранных моделей, которые вы используете, чтобы помочь выбрать лучшую модель.
Векторы E являются невязками, которые одного размера как данные об оценке.

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

Можно проверять, являются ли предполагаемые модели стабильными и обратимыми путем отображения свойства 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. Затем передайте различия в количестве предполагаемых параметров и loglikelihoods к 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
10³ ×

    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) любой:

Вызов forecast и передача последних нескольких строк YF
Симуляция нескольких временных рядов с simulate

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

Сгенерируйте 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);

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 и неуверенности по поводу направления уровней Казначейского векселя.

Документация