Моделирование путей выборки модели динамической регрессии с коммутацией Маркова
Y = simulate(Mdl,numObs,Name,Value)'NumPaths',1000,'Y0',Y0 моделирует 1000 sample paths и инициализирует динамический компонент каждой подмодели с помощью данных ответа предварительной выборки 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 ВВП 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.
Вычислить дистрибутив Monte Carlo
Для каждой переменной и пути вычислите среднее значение процесса.
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). (К. Милас, П. Ротман и Д. ван Дийк, ред.). Амстердам: Emerald Group Publishing Limited, 2006.
[2] Гамильтон, Дж. Д. «Анализ временных рядов, подверженных изменениям в режиме». Журнал эконометрики. Том 45, 1990, стр. 39-70.
[3] Гамильтон, Джеймс Д. Анализ временных рядов. Принстон, Нью-Джерси: Princeton University Press, 1994.