Этот пример показывает, как создать и оценить модель пространства состояний, содержащую изменяющиеся во времени параметры.
Предположим, что модель AR (2) и MA (1) содержат латентный процесс. Существует 50 периодов, и процесс MA (1) выпадает из модели на последние 25 периодов. Уравнение состояния для первых 25 периодов
Для последних 25 периодов уравнение состояния является
где и являются Гауссовыми со средним 0 и стандартным отклонением 1.
Сгенерируйте случайную серию из 50 наблюдений от и, принимая, что серия начинается с 1,5 и 1, соответственно.
T = 50; ARMdl = arima('AR',{0.7,-0.2},'Constant',0,'Variance',1); MAMdl = arima('MA',0.6,'Constant',0,'Variance',1); x0 = [1.5 1; 1.5 1]; rng(1); x = [simulate(ARMdl,T,'Y0',x0(:,1)),... [simulate(MAMdl,T/2,'Y0',x0(:,2));nan(T/2,1)]];
Последние 25 значений для моделируемых данных MA (1) отсутствуют.
Предположим далее, что латентные процессы измеряются с помощью
за первые 25 периодов, и
за последние 25 периодов. является Гауссовым со средним 0 и стандартным отклонением 1.
Сгенерируйте наблюдения, используя процесс случайного скрытого состояния (x
) и уравнение наблюдения.
y = 2*sum(x','omitnan')'+randn(T,1);
Вместе латентный процесс и уравнения наблюдений составляют модель пространства состояний. Предполагая, что коэффициенты являются неизвестными параметрами, модель пространства состояний является
Написание функции, которая задает, как параметры в params
сопоставить с матрицами модели пространства состояний, начальными значениями состояний и типом состояния.
% Copyright 2015 The MathWorks, Inc. function [A,B,C,D,Mean0,Cov0,StateType] = AR2MAParamMap(params,T) %AR2MAParamMap Time-variant state-space model parameter mapping function % % This function maps the vector params to the state-space matrices (A, B, % C, and D), the initial state value and the initial state variance (Mean0 % and Cov0), and the type of state (StateType). From periods 1 to T/2, the % state model is an AR(2) and an MA(1) model, and the observation model is % the sum of the two states. From periods T/2 + 1 to T, the state model is % just the AR(2) model. A1 = {[params(1) params(2) 0 0; 1 0 0 0; 0 0 0 params(3); 0 0 0 0]}; B1 = {[1 0; 0 0; 0 1; 0 1]}; C1 = {params(4)*[1 0 1 0]}; Mean0 = ones(4,1); Cov0 = 10*eye(4); StateType = [0 0 0 0]; A2 = {[params(1) params(2) 0 0; 1 0 0 0]}; B2 = {[1; 0]}; A3 = {[params(1) params(2); 1 0]}; B3 = {[1; 0]}; C3 = {params(5)*[1 0]}; A = [repmat(A1,T/2,1);A2;repmat(A3,(T-2)/2,1)]; B = [repmat(B1,T/2,1);B2;repmat(B3,(T-2)/2,1)]; C = [repmat(C1,T/2,1);repmat(C3,T/2,1)]; D = 1; end
Сохраните этот код в файле с именем AR2MAParamMap
и поместите его в путь MATLAB ®.
Создайте модель пространства состояний путем передачи функции AR2MAParamMap
как указатель на функцию, чтобы ssm
.
Mdl = ssm(@(params)AR2MAParamMap(params,T));
ssm
неявно задает модель пространства состояний. Обычно вы не можете проверить неявно определенные модели пространства состояний.
Передайте наблюдаемые отклики (y
) к estimate
для оценки параметров. Задайте положительные начальные значения для неизвестных параметров.
params0 = 0.1*ones(5,1); EstMdl = estimate(Mdl,y,params0)
Method: Maximum likelihood (fminunc) Sample size: 50 Logarithmic likelihood: -114.957 Akaike info criterion: 239.913 Bayesian info criterion: 249.473 | Coeff Std Err t Stat Prob --------------------------------------------------- c(1) | 0.47870 0.26634 1.79733 0.07229 c(2) | 0.00809 0.27179 0.02975 0.97626 c(3) | 0.55735 0.80958 0.68844 0.49118 c(4) | 1.62679 0.41622 3.90848 0.00009 c(5) | 1.90021 0.49563 3.83391 0.00013 | | Final State Std Dev t Stat Prob x(1) | -0.81229 0.46815 -1.73511 0.08272 x(2) | -0.31449 0.45918 -0.68490 0.49341 EstMdl = State-space model type: <a href="matlab: doc ssm">ssm</a> State vector length: Time-varying Observation vector length: 1 State disturbance vector length: Time-varying Observation innovation vector length: 1 Sample size supported by model: 50 State variables: x1, x2,... State disturbances: u1, u2,... Observation series: y1, y2,... Observation innovations: e1, e2,... State equations of period 1, 2, 3,..., 25: x1(t) = (0.48)x1(t-1) + (8.09e-03)x2(t-1) + u1(t) x2(t) = x1(t-1) x3(t) = (0.56)x4(t-1) + u2(t) x4(t) = u2(t) State equations of period 26: x1(t) = (0.48)x1(t-1) + (8.09e-03)x2(t-1) + u1(t) x2(t) = x1(t-1) State equations of period 27, 28, 29,..., 50: x1(t) = (0.48)x1(t-1) + (8.09e-03)x2(t-1) + u1(t) x2(t) = x1(t-1) Observation equation of period 1, 2, 3,..., 25: y1(t) = (1.63)x1(t) + (1.63)x3(t) + e1(t) Observation equation of period 26, 27, 28,..., 50: y1(t) = (1.90)x1(t) + e1(t) Initial state distribution: Initial state means x1 x2 x3 x4 1 1 1 1 Initial state covariance matrix x1 x2 x3 x4 x1 10 0 0 0 x2 0 10 0 0 x3 0 0 10 0 x4 0 0 0 10 State types x1 x2 x3 x4 Stationary Stationary Stationary Stationary
Предполагаемые параметры находятся в пределах 1 стандартной ошибки от их истинных значений, но стандартные ошибки довольно высоки. Поверхности вероятностей моделей пространства состояний могут содержать локальные максимумы. Поэтому рекомендуется попробовать несколько начальных значений параметров или рассмотреть возможность использования refine
.
estimate
| forecast
| refine
| simulate
| ssm