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

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

Загрузка и преобразование данных

Файл Data_USEconModel поставляется с программным обеспечением Econometrics Toolbox™. Файл содержит временные ряды из базы данных Федерального резервного банка данных экономики Сент-Луиса (FRED) в табличном массиве. Этот пример использует три временных рядов:

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

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

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

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 подходит для полных матриц для авторегрессивных матриц.

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

Для двух моделей 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

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

Чтобы найти лучшую модель в наборе, минимизируйте информационный критерий Акайке (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

Прогноз показывает увеличение реального ВВП и 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

Графики показывают рост ВВП, умеренный или незначительный рост M1 и неопределенность в отношении направления ставок T-векселей.