Создайте шаблон оценки модели в пространстве состояний

Чтобы оценить модель состояния-sapce, оценочная функция требует шаблона, задающего структуру модели, какие параметры являются неизвестными или фиксированными, и другими характеристиками.

Создайте модель в пространстве состояний для динамической факторной модели. Задайте значения NaN для неизвестных параметров и укажите, что состояния являются стационарными процессами путем установки их типа (StateType) к 0.

numStates = 2;
numSeries = 4;
A = tril(nan(numStates));
B = eye(numStates);
C = nan(numSeries,numStates);
D = diag(NaN(numSeries,1));
Mdl = ssm(A,B,C,D,StateType=[0 0]);

Mdl ssm объект, представляющий шаблон для оценки.

Загрузите и предварительно обработайте данные

Загрузите США макроэкономический набор данных Data_USEconModel.mat. Данные содержат ежеквартальные измерения 14 экономических рядов от Q1:1947 до Q1:2009.

load Data_USEconModel

DataTable 249 14 расписание, которое содержит данные и демонстрационные метки времени. Для получения дополнительной информации введите Description в командной строке.

Извлеките CPI, GDP, 3-месячный доход по казначейским векселям и уровень безработицы.

DataLevel = DataTable{:,["CPIAUCSL" "GDP" "TB3MS" "UNRATE"]};
anyMissing = any(any(ismissing(DataLevel)))

anyMissing = logical
   1

DataLevel 249 4 числовая матрица измерений. По крайней мере один ряд содержит отсутствующие значения, но Фильтр Калмана вмещает недостающие наблюдения путем продвижения вперед распределения текущего состояния как следующего прогноза состояния.

Преобразуйте уровни в возвраты путем передачи данных price2ret.

DataGrowth = price2ret(DataLevel)*100;

Поскольку преобразование использует первое различие, DataGrowth имеет тот меньше наблюдения, чем DataLevel.

Зарезервируйте итоговые 3 года (12 четвертей) данных как выборка затяжки. Извлеките соответствующие даты.

numPeriods = 12;
Y = DataGrowth(1:(end-numPeriods),:);       % In-sample data 
Yf = DataGrowth((end-numPeriods+1):end,:);  % Holdout sample
T = size(DataGrowth,1);

dates = DataTable.Time(2:end-numPeriods);
datesf = DataTable.Time(end-numPeriods+1:end);

Оценочная модель

Процедура оценки требует начальных значений для всех неизвестных параметров. Создайте вектор из содержания начальных значений. Например, если A0, C0, и D0 матрицы, содержащие начальные значения для AB, и C, расположите начальные значения путем создания векторного [A0(:,1); A0(:,2); C0(:,1); C0(:,2); D0(:,1); D0(:,2); D0(:,3); D0(:,4)]. Только неизвестные параметры требуют начальных значений.

params0 = [0.8; 0; 0.5; ones(numStates*numSeries,1); 0.5*std(Y,"omitNan")'];

Подбирайте модель к в выборочных данных. Задайте начальные значения и одномерную обработку многомерной модели. Выключите сводные данные оценки.

EstMdl = estimate(Mdl,Y,params0,Univariate=1,Display="off")

EstMdl = 
State-space model type: ssm

State vector length: 2
Observation vector length: 4
State disturbance vector length: 2
Observation innovation vector length: 4
Sample size supported by model: Unlimited

State variables: x1, x2,...
State disturbances: u1, u2,...
Observation series: y1, y2,...
Observation innovations: e1, e2,...

State equations:
x1(t) = (0.60)x1(t-1) + u1(t)
x2(t) = -(0.23)x1(t-1) + (0.99)x2(t-1) + u2(t)

Observation equations:
y1(t) = -(0.07)x1(t) + (0.12)x2(t) + (0.61)e1(t)
y2(t) = -(0.70)x1(t) + (0.18)x2(t) + (0.63)e2(t)
y3(t) = -(6.87)x1(t) + (0.07)x2(t) + (16.18)e3(t)
y4(t) = (5.30)x1(t) + (0.07)x2(t) + (3.96)e4(t)

Initial state distribution:

Initial state means
 x1  x2 
  0   0 

Initial state covariance matrix
     x1     x2    
 x1   1.57  -0.53 
 x2  -0.53  76.34 

State types
     x1          x2     
 Stationary  Stationary 

Проверяйте, устойчива ли предполагаемая модель состояния VAR.

varlagop = LagOp({eye(numStates) EstMdl.A});
isStable(varlagop)

ans = logical
   1

Модель состояния VAR устойчива.

Симулируйте асинхронный релиз измерений

Чтобы симулировать асинхронный релиз экономических измерений в четверти, предположите, что каждое измерение выпущено в последовательные недели. Для наблюдения t в выборке затяжки:

Сложите каждую матрицу.

YfStretch = nan(numPeriods*numSeries,numSeries);
for t = 1:numPeriods
    for j = 1:numSeries
        YfStretch((t-1)*numSeries + j,j) = Yf(t,j);
    end
end

YfStretch numPeriods*numSeries- numSeries (48 4) блочная матрица асинхронно выпущенных измерений за период расчета затяжки. Например, блоком 1 является YfStretch(1:4,:) и это представляет первое наблюдение в выборке затяжки Yf(1,:) где эти четыре измерения выпущены в разное время или подпериоды (строки), блоком 2 является YfStretch(5:8,:) и это представляет второе наблюдение в выборке затяжки Yf(2,:) где четыре измерения выпущены в разное время и так далее.

Предполагаемая модель EstMdl калибруется так, чтобы переход состояний каждый период (строка в данных). Состояния в выборке затяжки должны перейти только после обработки входящего темпа роста безработицы. Чтобы достигнуть этого действия, создайте модель для наблюдений подпериода так, чтобы состояния остались то же самое, когда это обрабатывает наблюдение. Символически,

с вырожденным распределением начального состояния.

EstMdlAux = ssm(eye(numStates),zeros(numStates),EstMdl.C,EstMdl.D,...
    Mean0=zeros(numStates,1),Cov0=zeros(numStates));

Nowcast модель

Получите моменты распределенности в конце в периоде расчета.

[states,~,output] = filter(EstMdl,Y);
currentState = states(end,:)';
CurrentStateCov = output(end).FilteredStatesCov;

Предварительно выделите массивы для nowcasts и даты расширенных периодов.

newStateF = nan(numPeriods*numSeries,numStates);
NewStateCovF = cell(numPeriods*numSeries,1);
datesfStretch = repmat(dates(end),numPeriods*numSeries,1);

Теперь бросьте модель при помощи этой процедуры для каждой блочной матрицы в расширенной выборке затяжки:

tStretch = 0;   % Index of stretched forecast horizon
wk = 0;         

for t = 1:numPeriods
    for j = 1:numSeries
        tStretch = tStretch + 1;

        if j == 1   % Advance to next full period
            [newStateF(tStretch,:),NewStateCovF{tStretch}] = update(EstMdl,YfStretch(tStretch,:),...
                currentState,CurrentStateCov);
            wk = 0;
            datesfStretch(tStretch) = datesf(t);

        else        % Stretch current period
            [newStateF(tStretch,:),NewStateCovF{tStretch}] = update(EstMdlAux,YfStretch(tStretch,:),...
                currentState,CurrentStateCov);
            wk = wk + 1;
            datesfStretch(tStretch) = datesf(t) + wk;
        end

        currentState = newStateF(tStretch,:)';
        CurrentStateCov = NewStateCovF{tStretch};
    end
end

newStateF 48 2 блочная матрица требований nowcasted и предоставлений в подпериоды периода затяжки. NewStateCovF 48 1 массив ячеек ковариационных матриц 2 на 2 для каждого nowcast.

Постройте оценки спроса и предложения в выборке в течение периодов вне 1 993 и постройте nowcasts.

dtp = dates > datetime(1994,1,1);

figure
h = plot(dates(dtp),states(dtp,:),"-",datesfStretch,newStateF,"--");
hold on
h(3).Color = h(1).Color;
h(4).Color = h(2).Color;
yl = ylim;
hp = patch([datesfStretch(1) datesfStretch(end) datesfStretch(end) datesfStretch(1)],...
    yl([1,1,2,2]),[0.8 0.8 0.8]);
uistack(hp,'bottom');
legend([h; hp],...
    ["Demand (in-sample)"; "Supply (in-sample)"; "Demand (updated)"; "Supply (updated)"; "Forecast horizon"],...
    Location="best")
hold off