Симулируйте демонстрационные пути переключающей Маркова модели динамической регрессии
Y = simulate(Mdl,numObs,Name,Value)'NumPaths',1000,'Y0',Y0 симулирует 1000 демонстрационные пути и инициализируют динамический компонент каждой подмодели при помощи преддемонстрационных данных об ответе Y0.
[ также возвращает симулированные инновационные пути Y,E,StatePaths] = simulate(___)E и симулированные пути состояния StatePaths, использование любой из комбинаций входных аргументов в предыдущих синтаксисах.
Симулируйте путь к ответу из переключающей Маркова модели динамической регрессии с двумя состояниями для 1D процесса ответа. Этот пример использует произвольные значения параметров.
Создайте полностью заданную модель
Создайте модель дискретной цепи Маркова с двумя состояниями, которая описывает механизм переключения режима. Пометьте режимы.
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")
Полагайте, что модель в Симулирует Путь к Ответу.
Создайте переключающую Маркова модель динамической регрессии.
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
YE, и 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)
Считайте переключающую Маркова модель динамической регрессии с двумя состояниями послевоенных США действительным темпом роста GDP. Модели представили оценки параметра в [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];
Создайте переключающую Маркова модель динамической регрессии, которая описывает поведение темпа роста GDP США.
Mdl = msVAR(mc,mdl);
Mdl полностью заданный msVAR объект.
Сгенерируйте один случайный путь 100 ответов, соответствующих инноваций и состояний из модели.
rng(1) % For reproducibility
[y,e,sp] = simulate(Mdl,100);y 100 1 вектор из уровней GDP и e 100 1 вектор из соответствующих инноваций. sp 100 1 вектор из индексов состояния.
Полагайте, что переключающая Маркова модель в Симулирует Уровни GDP США и Экономические состояния, но принимает, что подмодели являются 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
Сгенерируйте случайные пути из переключающей Маркова модели динамической регрессии с тремя состояниями для 2D процесса ответа 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);Симулируйте разнообразные пути
Сгенерируйте 1 000 случайных путей ответов для 50 временных шагов. Запустите все симуляции в первом состоянии.
rng(10); % For reproducibility S0 = [1 0 0]; Y = simulate(Mdl,50,'S0',S0,'NumPaths',1000);
Y 50 2 1 000 массивов путей к симулированному отклику.
Вычислите распределение Монте-Карло
Для каждой переменной и пути, вычислите среднее значение процесса.
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] Chauvet, M. и Дж. Д. Гамильтон. "Датируя Поворотные моменты Делового цикла". В Нелинейном Анализе Деловых циклов (Вклады в Экономический анализ, Объем 276). C. Милас, П. Ротмен, и Д. ван Дейк, редакторы). Амстердам: Emerald Group Publishing Limited, 2006.
[2] Гамильтон, J. D. "Анализ Временных рядов Согласно Изменениям в Режиме". Журнал Эконометрики. Издание 45, 1990, стр 39–70.
[3] Гамильтон, анализ временных рядов Джеймса Д. Принстон, NJ: Издательство Принстонского университета, 1994.