function [A,B,C,D,mean0,cov0,stateType,deflatedData] = Example_YieldsExogenous(param,Yields,Macro,maturities)
%EXAMPLE_YIELDSEXOGENOUS Parameter-to-matrix map for Diebold-Li model
%augmented with exogenous macroeconomic factors
%
% Syntax:
%
%   [A,B,C,D,mean0,cov0,stateType,deflatedYield] = ...
%                                 Example_YieldsExogenous(param,Yields,Macro,maturities)
%
% Note:
%
%   To create an ssm model object that implicitly maps Diebold-Li model
%   parameters to SSM coefficients, parameterize the function as an
%   anonymous function that accepts a single column vector of estimated
%   parameters:
%
%   f = @(param)Example_YieldsExogenous(param,Yields,Macro,maturities)
%
%   Example_YieldsMacro supports the Diebold-Li example Demo_DieboldLiModel;
%   the function can change in a future release without notice.
%
% Description:
%
%   For observation vector y(t) and state vector x(t), create a state-space 
%   model (SSM) of the form:
%
%   State equation:       x(t) = A * x(t-1) + B * u(t)
%   Observation equation: y(t) = C * x(t)   + D * e(t)
%
%   where u(t) and e(t) are uncorrelated, unit-variance white noise vector
%   processes.  In the example, the 6-by-1 state vector x(t) represents the
%   three Diebold-Li model factors and three macroeconomic factors, and A
%   is the 6-by-6 state transition matrix.  The dimension of the observed
%   yield vector, y(t), is determined by the number of maturities included
%   in the time series of yield curves.
%
%   This function imposes the following constraints:
%   * The macroeconomic variables are exogenous (A(4:6,1:3) = zeros(3)).
%   * The covariance matrix of u(t), Q = B*B', is diagonal.
%   * The covariance matrix of e(t), H = D*D', is diagonal.  
%   This function also accounts for the means of the six factors (states,
%   x(t)) and deflates the input observed yields y(t) to estimate the
%   factor offsets.
%

% Copyright 2021 The MathWorks, Inc.
    numBonds = size(Yields,2);
    numMacro = size(Macro,2);
    numStates = 3 + numMacro;

    lambda = param(end);
    lambdaTau = lambda.*maturities;
    expLambdaTau = exp(-lambdaTau);
    loading = [ones(numBonds,1), (1 - expLambdaTau)./lambdaTau, ...
        (1 - expLambdaTau)./lambdaTau - expLambdaTau];
    C = blkdiag(loading,eye(numMacro));

    mask = true(numStates);
    mask(4:end,1:3) = false;
    A = zeros(numStates); 
    A(mask) = param(1:sum(mask(:)));    

    iOffset = sum(mask(:));
    B = diag(param((iOffset + 1):(iOffset + numStates)));

    iOffset = iOffset + numStates;
    stdev = param((iOffset + 1):(iOffset + numBonds));
    D = [diag(stdev);zeros(numMacro,numBonds)];

    iOffset = iOffset + numBonds;
    mu = param((iOffset + 1):(iOffset + numStates));
    deflatedData = [Yields,Macro] - (C*mu)';

    mean0     = [];
    cov0      = [];
    stateType = [];
end
