В этом примере показан анализ модели VAR.
В этом разделе приведен пример рабочего процесса, описанного в разделе Рабочий процесс модели VAR. В примере используются три временных ряда: ВВП, M1 денежная масса и 3-месячная ставка Т-векселя. В примере показано:
Загрузка и преобразование данных для стационарности
Разделение преобразованных данных на интервалы предварительной выборки, оценки и прогноза для поддержки эксперимента по тестированию
Создание нескольких моделей
Подбор моделей к данным
Решение о том, какая из моделей является лучшей
Составление прогнозов на основе наилучшей модели
Файл Data_USEconModel поставляется с программным обеспечением Econometrics Toolbox™. Файл содержит временные ряды из базы данных Federal Reserve Bank of St. Louis Economics Data (FRED) в табличном массиве. В этом примере используются три временных ряда:
ВВП (GDP)
M1 денежная масса (M1SL)
3-месячная ставка Т-счета (TB3MS)
Загрузите набор данных. Создайте переменную для реального ВВП.
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-bill не показывают экспоненциального роста. Чтобы противостоять тенденциям в реальном ВВП и M1, возьмем разницу логарифмов данных. Кроме того, стабилизируйте серию T-bill, взяв первую разницу. Синхронизируйте ряды дат таким образом, чтобы данные имели одинаковое количество строк для каждого столбца.
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
Для данных можно выбрать множество различных моделей. В этом примере используются четыре модели.
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 подгоняет полные матрицы для авторегрессивных матриц.
Чтобы оценить качество моделей, создайте векторы индексов, которые разделяют данные ответа на три периода: предварительную выборку, оценку и прогноз. Поместите модели в оценочные данные, используя период предварительной выборки для предоставления запаздывающих данных. Сравните прогнозы подходящих моделей с данными прогноза. Период оценки находится в выборке, а период прогноза находится вне выборки (также известный как обратное тестирование).
Для двух моделей 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 значения являются логическими связями подогнанных моделей, которые используются для выбора наилучшей модели.
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. Затем передайте различия в количестве оцениваемых параметров и логарифмов в 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) является лучшей моделью.
Чтобы найти лучшую модель в наборе, минимизируйте информационный критерий Акайке (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
Прогноз показывает рост реального ВВП и 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, а также неопределенность в отношении направления ставок Т-векселя.
aicbic | estimate | forecast | lratiotest | simulate