Этот пример показывает, как анализировать модель 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
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) любой:
Вызов 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);
Преобразуйте прогнозы путем отмены масштабирования, и дифференцирование применилось к исходным данным. Убедитесь, что вставили последнее наблюдение в начале временных рядов перед использованием 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 и неуверенности по поводу направления уровней Казначейского векселя.
aicbic | estimate | forecast | lratiotest | simulate