Симулируйте пути расчета модели динамической регрессии переключения Маркова
Y = simulate(Mdl,numObs,Name,Value)'NumPaths',1000,'Y0',Y0 моделирует 1000 выборочные пути и инициализирует динамический компонент каждой подмодели с помощью предварительных примеров данных отклика Y0.
[ также возвращает моделируемые пути инноваций Y,E,StatePaths] = simulate(___)E и моделируемые пути состояния StatePaths, использование любой комбинации входных аргументов в предыдущих синтаксисах.
Симулируйте путь отклика из модели динамической регрессии с двумя состояниями Маркова-переключения для процесса 1-D отклика. Этот пример использует произвольные значения параметров.
Создайте полностью заданную модель
Создайте модель дискретной цепи Маркова с двумя состояниями, которая описывает механизм переключения режима. Пометьте режимы.
P = [0.9 0.1; 0.3 0.7]; mc = dtmc(P,'StateNames',["Expansion" "Recession"]);
mc является полностью заданным dtmc объект.
Для каждого режима используйте arima для создания модели AR, которая описывает процесс ответа в режиме. Сохраните подмодели в векторе.
mdl1 = arima('Constant',5,'AR',[0.3 0.2],... 'Variance',2); mdl2 = arima('Constant',-5,'AR',0.1,... 'Variance',1); mdl = [mdl1; mdl2];
mdl1 и mdl2 полностью заданы arima объекты.
Создайте модель динамической регрессии переключения Маркова из механизма переключения mc и вектор подмоделей mdl.
Mdl = msVAR(mc,mdl);
Mdl является полностью заданным msVAR объект.
Моделируйте путь отклика
Сгенерируйте один путь случайного отклика длины 50 из модели.
rng(1); % For reproducibility
y = simulate(Mdl,50);y является вектором 50 на 1 из одного пути отклика.
Постройте график пути ответа.
figure plot(y) xlabel("Time") ylabel("Response")
Рассмотрим модель в Simulate Response Path.
Создайте модель динамической регрессии с переключением Маркова.
P = [0.9 0.1; 0.3 0.7]; mc = dtmc(P,'StateNames',["Expansion" "Recession"]); mdl1 = arima('Constant',5,'AR',[0.3 0.2],... 'Variance',2); mdl2 = arima('Constant',-5,'AR',0.1,... 'Variance',1); mdl = [mdl1; mdl2]; Mdl = msVAR(mc,mdl);
Симулируйте 3 пути отклика, инноваций и индекса состояний из 5 наблюдений от модели.
rng('default') % For reproducibility [Y,E,SP] = simulate(Mdl,5,'NumPaths',3)
Y = 5×3
-5.7605 10.9496 11.4633
-5.7002 8.5772 -3.1268
4.2446 10.7774 -5.6161
-3.1665 -2.2920 -5.2677
-3.8995 -4.7403 -6.3141
E = 5×3
-0.2050 0.9496 1.4633
-0.1241 -1.7076 0.7269
2.1068 1.0143 -0.3034
1.4090 1.6302 0.2939
1.4172 0.4889 -0.7873
SP = 5×3
2 1 1
2 1 2
1 1 2
2 2 2
2 2 2
Y, E, и SP представляют собой матрицы 5 на 3. Столбцы представляют отдельные, независимые пути.
Симулируйте один путь откликов, инноваций и состояний в горизонт симуляции длины 50. Затем постройте график каждого пути отдельно.
[y,e,sp] = simulate(Mdl,50); figure subplot(3,1,1) plot(y) ylabel('Response') grid on subplot(3,1,2) plot(e) ylabel('Innovation') grid on subplot(3,1,3) plot(sp,'m') ylabel('State') yticks([1 2]) yticklabels(Mdl.StateNames)
Рассмотрим двухгосударственную марковскую динамическую регрессионую модель послевоенного реального темпа роста ВВП США. Модель имеет оценки параметров, представленные в [1].
Создайте модель дискретной цепи Маркова, которая описывает механизм переключения режима. Пометьте режимы.
P = [0.92 0.08; 0.26 0.74]; mc = dtmc(P,'StateNames',["Expansion" "Recession"]);
Создайте отдельные модели AR (0) (только постоянные) для этих двух режимов.
sigma = 3.34; % Homoscedastic models across states mdl1 = arima('Constant',4.62,'Variance',sigma^2); mdl2 = arima('Constant',-0.48,'Variance',sigma^2); mdl = [mdl1 mdl2];
Создайте модель динамической регрессии переключения Маркова, которая описывает поведение темпов роста ВВП США.
Mdl = msVAR(mc,mdl);
Mdl является полностью заданным msVAR объект.
Сгенерируйте один случайный путь из 100 откликов, соответствующих инноваций и состояний из модели.
rng(1) % For reproducibility
[y,e,sp] = simulate(Mdl,100);y является вектором скорости ВВП 100 на 1, и e является вектором 100 на 1 соответствующих инноваций. sp является вектором индексов состояний 100 на 1.
Рассмотрим модель марковского переключения в Simulate US GDP Rates and Economic States, но предположим, что подмодели являются AR (1) вместо этого. Рассмотрим подбор кривой модель к наблюдениям в период 1960: Q1-2004: Q2.
Создайте шаблон модели для оценки. Укажите подмодели AR (1).
mc = dtmc(NaN(2),'StateNames',["Expansion" "Recession"]); ar1 = arima(1,0,0); Mdl = msVAR(mc,[ar1; ar1]);
Поскольку подмодели являются AR (1), каждому требуется одно предварительное наблюдение, чтобы инициализировать его динамический компонент для оценки .
Создайте модель, содержащую начальные значения параметров для процедуры оценки.
mc0 = dtmc(0.5*ones(2),'StateNames',["Expansion" "Recession"]); submdl01 = arima('Constant',1,'Variance',1,'AR',0.001); submdl02 = arima('Constant',-1,'Variance',1,'AR',0.001); Mdl0 = msVAR(mc0,[submdl01; submdl02]);
Загрузите данные. Преобразуйте весь набор в годовую последовательность скоростей.
load Data_GDP
qrate = diff(Data)./Data(1:(end - 1));
arate = 100*((1 + qrate).^4 - 1);Идентифицируйте предварительную выборку и оценку периодов дискретизации с помощью дат, связанных с годовым рядом ставок. Поскольку преобразование применяет первое различие, необходимо удалить первую дату наблюдения из исходной выборки.
dates = datetime(dates(2:end),'ConvertFrom','datenum',... 'Format','yyyy:QQQ','Locale','en_US'); estPrd = datetime(["1960:Q2" "2004:Q2"],'InputFormat','yyyy:QQQ',... 'Format','yyyy:QQQ','Locale','en_US'); idxEst = isbetween(dates,estPrd(1),estPrd(2)); idxPre = dates < estPrd(1);
Подбор модели к оценочным выборочным данным. Задайте предварительный образец наблюдения.
y0 = arate(idxPre);
EstMdl = estimate(Mdl,Mdl0,arate(idxEst),'Y0',y0);Симулируйте путь отклика от подобранной модели в течение периода оценки. Задайте предварительный образец наблюдения.
rng(1) % For reproducibility numObs = sum(idxEst); aratesim = simulate(EstMdl,numObs,'Y0',y0);
Постройте графики наблюдений и моделируемого пути годовых скоростей и идентифицируйте периоды рецессии с помощью recessionplot.
figure; plot(dates(idxEst),[arate(idxEst) aratesim]) recessionplot xlabel("Time") ylabel("Annualized GDP Rate") legend("Observed","Simulated");
Оцените точность оценки, используя моделируемые данные из известного процесса генерации данных (DGP). Этот пример использует произвольные значения параметров.
Создайте модель для DGP
Создайте полностью заданную модель дискретной цепи Маркова двух состояний для механизма переключения.
P = [0.7 0.3; 0.1 0.9]; mc = dtmc(P);
Для каждого состояния создайте полностью заданную модель AR (1) для процесса отклика.
% Constants C1 = 5; C2 = -2; % Autoregression coefficients AR1 = 0.4; AR2 = 0.2; % Variances V1 = 4; V2 = 2; % AR Submodels dgp1 = arima('Constant',C1,'AR',AR1,'Variance',V1); dgp2 = arima('Constant',C2,'AR',AR2,'Variance',V2);
Создайте полностью заданную модель динамической регрессии Маркова-переключения для DGP.
DGP = msVAR(mc,[dgp1,dgp2]);
Моделирование путей отклика от DGP
Сгенерируйте 10 путей случайного отклика длиной 1000 от DGP.
rng(1); % For reproducibility N = 10; n = 1000; Data = simulate(DGP,n,'Numpaths',N);
Data является 1000 на 10 матрицей симулированных откликов.
Создайте модель для оценки
Создайте частично заданную модель динамической регрессии Маркова-переключения, которая имеет ту же структуру, что и процесс генерации данных, но задает неизвестную матрицу перехода и неизвестные коэффициенты подкоэффициентов модели.
PEst = NaN(2); mcEst = dtmc(PEst); mdl = arima(1,0,0); Mdl = msVAR(mcEst,[mdl; mdl]);
Создайте модель, содержащую начальные значения
Создайте полностью заданную модель динамической регрессии Маркова-переключения, которая имеет ту же структуру, что и Mdl, но установите все оцениваемые параметры в начальные значения.
P0 = 0.5*ones(2); mc0 = dtmc(P0); mdl01 = arima('Constant',1,'AR',0.5,'Variance',2); mdl02 = arima('Constant',-1,'AR',0.5,'Variance',1); Mdl0 = msVAR(mc0,[mdl01,mdl02]);
Оценка моделей
Подбор модели к каждому моделируемому пути. Для каждого пути постройте график логарифмической правдоподобности при каждой итерации алгоритма EM.
c1 = zeros(N,1); c2 = zeros(N,1); v1 = zeros(N,1); v2 = zeros(N,1); ar1 = zeros(N,1); ar2 = zeros(N,1); PStack = zeros(2,2,N); figure hold on for i = 1:N EstModel = estimate(Mdl,Mdl0,Data(:,i),'IterationPlot',true); c1(i) = EstModel.Submodels(1).Constant; c2(i) = EstModel.Submodels(2).Constant; v1(i) = EstModel.Submodels(1).Covariance; v2(i) = EstModel.Submodels(2).Covariance; ar1(i) = EstModel.Submodels(1).AR{1}; ar2(i) = EstModel.Submodels(2).AR{1}; PStack(:,:,i) = EstModel.Switch.P; end hold off
Оцените точность
Вычислите среднее значение Монте-Карло для каждого предполагаемого параметра.
c1Mean = mean(c1); c2Mean = mean(c2); v1Mean = mean(v1); v2Mean = mean(v2); ar1Mean = mean(ar1); ar2Mean = mean(ar2); PMean = mean(PStack,3);
Сравните параметры населения с соответствующими оценками Монте-Карло.
DGPvsEstimate = [...
C1 c1Mean
C2 c2Mean
V1 v1Mean
V2 v2Mean
AR1 ar1Mean
AR2 ar2Mean]DGPvsEstimate = 6×2
5.0000 5.0260
-2.0000 -1.9615
4.0000 3.9710
2.0000 1.9903
0.4000 0.4061
0.2000 0.2017
P
P = 2×2
0.7000 0.3000
0.1000 0.9000
PEstimate = PMean
PEstimate = 2×2
0.7065 0.2935
0.1023 0.8977
Сгенерируйте случайные пути из модели динамической регрессии Маркова с тремя состояниями для 2-D процесса отклика VARX. Этот пример использует произвольные значения параметров для DGP.
Создайте полностью заданную модель для DGP
Создайте модель дискретной цепи Маркова трех состояний для механизма переключения.
P = [10 1 1; 1 10 1; 1 1 10]; mc = dtmc(P);
mc является полностью заданным dtmc объект. dtmc нормализует строки P чтобы они суммировали 1.
Для каждого режима используйте varm для создания модели VAR, которая описывает процесс отклика в режиме. Задайте все значения параметров.
% Constants C1 = [1;-1]; C2 = [2;-2]; C3 = [3;-3]; % Autoregression coefficients AR1 = {}; AR2 = {[0.5 0.1; 0.5 0.5]}; AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]}; % Regression coefficients Beta1 = [1;-1]; Beta2 = [2 2;-2 -2]; Beta3 = [3 3 3;-3 -3 -3]; % Innovations covariances Sigma1 = [1 -0.1; -0.1 1]; Sigma2 = [2 -0.2; -0.2 2]; Sigma3 = [3 -0.3; -0.3 3]; % VARX submodels mdl1 = varm('Constant',C1,'AR',AR1,'Beta',Beta1,'Covariance',Sigma1); mdl2 = varm('Constant',C2,'AR',AR2,'Beta',Beta2,'Covariance',Sigma2); mdl3 = varm('Constant',C3,'AR',AR3,'Beta',Beta3,'Covariance',Sigma3); mdl = [mdl1; mdl2; mdl3];
mdl содержит три полностью заданных varm объекты модели.
Для DGP создайте полностью заданную модель динамической регрессии Маркова-переключения из механизма переключения mc и подмодели mdl.
Mdl = msVAR(mc,mdl);
Mdl является полностью заданным msVAR модель.
Симулируйте данные, игнорирующие регрессионый компонент
Если вы не поставляете экзогенные данные, simulate игнорирует компоненты регрессии в подмоделях. Симулируйте 3 отдельных, независимых пути откликов, инноваций и индексов состояний длины 5 от модели.
rng(1); % For reproducibility [Y,E,SP] = simulate(Mdl,5,'NumPaths',3)
Y =
Y(:,:,1) =
5.2387 1.5297
4.4290 4.2738
1.1668 -1.2905
-0.9654 -0.2028
-0.2701 0.8993
Y(:,:,2) =
2.7737 -2.5383
-0.8651 -1.1046
-0.0511 0.3696
0.5826 -0.8926
2.4022 -0.6912
Y(:,:,3) =
3.5443 0.8768
4.9748 -0.7956
5.7213 0.8073
4.2473 0.5805
2.7972 -1.3340
E =
E(:,:,1) =
1.2387 1.5297
-0.3434 2.8896
0.1668 -0.2905
-1.9654 0.7972
-1.2701 1.8993
E(:,:,2) =
1.7737 -1.5383
-1.8651 -0.1046
-1.0511 1.3696
-0.4174 0.1074
1.4022 0.3088
E(:,:,3) =
-0.4557 0.8768
1.1150 -1.0061
1.3134 0.7176
-0.6941 -0.6838
1.7972 -0.3340
SP = 5×3
2 1 2
2 1 2
1 1 2
1 1 2
1 1 1
Y и E представляют собой массивы симулированных откликов и инноваций 5 на 2 на 3 соответственно. Строки соответствуют временным точкам, столбцы соответствуют переменным в системе, а страницы соответствуют путям. SP - матрица 5 на 3, столбцы которой соответствуют путям.
Симулируйте один путь откликов, инноваций и состояний в горизонт симуляции длины 50. Затем постройте график каждого пути отдельно.
rng0 = rng; % Store settings to reproduce state sequence. [Y,E,SP] = simulate(Mdl,50); figure subplot(3,1,1) plot(Y) ylabel("Response") grid on legend(["y_1" "y_2"]) subplot(3,1,2) plot(E) ylabel("Innovation") grid on legend(["e_1" "e_2"]) subplot(3,1,3) plot(SP,'m') ylabel("State") yticks([1 2 3])
Моделируйте данные, включая регрессионный компонент
Симулируйте экзогенные данные для трех регрессоров, сгенерировав 50 случайных наблюдений из 3-D стандартного Гауссова распределения.
X = randn(50,3);
Сгенерируйте один случайный ответ, инновации и путь состояния длины 50. Задайте моделируемые экзогенные данные для компонентов регрессии подмодели. Постройте график результатов.
rng(rng0); % Reproduce state sequence in previous simulation. [Y,E,SP] = simulate(Mdl,50,'X',X); figure subplot(3,1,1) plot(Y) ylabel("Response") grid on legend(["y_1" "y_2"]) subplot(3,1,2) plot(E) ylabel("Innovation") grid on legend(["e_1" "e_2"]) subplot(3,1,3) plot(SP,'m') ylabel("State") yticks([1 2 3])
Рассмотрите модель в Симуляции путей из Модели с Подмоделями VARX.
Создайте полностью заданную модель
Создайте модель переключения Маркова, исключая регрессионый компонент.
P = [10 1 1; 1 10 1; 1 1 10];
mc = dtmc(P);
C1 = [1;-1];
C2 = [2;-2];
C3 = [3;-3];
AR1 = {};
AR2 = {[0.5 0.1; 0.5 0.5]};
AR3 = {[0.25 0; 0 0] [0 0; 0.25 0]};
Sigma1 = [1 -0.1; -0.1 1];
Sigma2 = [2 -0.2; -0.2 2];
Sigma3 = [3 -0.3; -0.3 3];
mdl1 = varm('Constant',C1,'AR',AR1,'Covariance',Sigma1);
mdl2 = varm('Constant',C2,'AR',AR2,'Covariance',Sigma2);
mdl3 = varm('Constant',C3,'AR',AR3,'Covariance',Sigma3);
mdl = [mdl1; mdl2; mdl3];
Mdl = msVAR(mc,mdl);Моделирование нескольких путей
Сгенерируйте 1000 случайных путей откликов для 50 временных шагов. Запустите все симуляции в первом состоянии.
rng(10); % For reproducibility S0 = [1 0 0]; Y = simulate(Mdl,50,'S0',S0,'NumPaths',1000);
Y представляет собой массив 50 на 2 на 1000 симулированные отклики путей.
Вычисление распределения Монте-Карло
Для каждой переменной и пути вычислите среднее значение процесса.
mus = mean(Y,1);
Для каждой переменной постройте график распределения Монте-Карло среднего процесса.
figure h1 = histogram(mus(1,1,:),'Normalization',"probability",... 'BinWidth',0.1); hold on h2 = histogram(mus(1,2,:),'Normalization',"probability",... 'BinWidth',0.1); legend(["y_1" "y_2"]) title('Process Means') hold off
Для каждой переменной и пути вычислите стандартное отклонение процесса.
sigmas = std(Y,0,1);
Для каждой переменной постройте график распределения Монте-Карло стандартного отклонения процесса.
figure h1 = histogram(sigmas(1,1,:),'Normalization',"probability",... 'BinWidth',0.05); hold on h2 = histogram(sigmas(1,2,:),'Normalization',"probability",... 'BinWidth',0.05); legend(["y_1" "y_2"]) title('Process Standard Deviations') hold off
[1] Шове, М. и Дж. Д. Гамильтон. «Поворотные точки делового цикла знакомств». В нелинейном анализе бизнес-циклов (вклад в экономический анализ, том 276). (C. Milas, P. Rothman, and D. van Dijk, eds.). Амстердам: Emerald Group Publishing Limited, 2006.
[2] Гамильтон, Дж. Д. «Анализ временных рядов, подверженных изменениям в режиме». Журнал эконометрики. Том 45, 1990, стр. 39-70.
[3] Гамильтон, Джеймс Д. Анализ временных рядов. Princeton, NJ: Princeton University Press, 1994.