%CESLDUR.M
%  This program solves the expenditure share problem discussed in
%  Section 6 of Binder and Pesaran (1999) using the LDU factorization
%  procedure described in Proposition 5.1 of Binder and Pesaran (1999). 
%
%  The program returns the optimal expenditure shares for T-t+1 periods,
%  where the first period, t, and the terminal period, T, can be chosen 
%  by the user.
%  There are three commodity groups in the model considered in this
%  program, and the adjustment costs are of first order. The program
%  is based on the canonical form (6.14) in Binder and Pesaran (1999). 
%
%  If you find any error in this program, please send e-mail to: 
%  binder@glue.umd.edu
%  Copyright: Michael Binder and M. Hashem Pesaran
%  Current Version: 03/02/2000
%
%  The paper 'Solution of Finite-Horizon Multivariate Linear 
%  Rational Expectations Models and Sparse Linear Systems' 
%  describing the theory underlying this program can be 
%  downloaded from:
%  http://www.inform.umd.edu/econ/mbinder
%  and
%  http://www.econ.cam.ac.uk/faculty/pesaran
%
%  See the user notes at the first of these two URL's before using 
%  this program.

clear variables 

% Specify Parameter Values (Notation as in Binder and Pesaran, 1999)
beta = .96;
alpha = zeros(3,1);
alpha(1,1) = .2994;
alpha(2,1) = .5053;
alpha(3,1) = 1-alpha(1,1)-alpha(2,1);
delta = zeros(3,1);
delta(1,1) = -.2361;
delta(2,1) = .3228;
delta(3,1) = -delta(1,1)-delta(2,1);
omega0 = zeros(3,1);
omega0(1,1) = .5384;
omega0(2,1) = .3054;
omega0(3,1) = .1562;
Gamma = zeros(3,3);
Gamma(1,1) = .2695;
Gamma(1,2) = -.1154;
Gamma(1,3) = -.1541;
Gamma(2,1) = Gamma(1,2);
Gamma(2,2) = -.0420;
Gamma(2,3) = .1575;
Gamma(3,1) = Gamma(1,3);
Gamma(3,2) = Gamma(2,3);
Gamma(3,3) = -Gamma(3,1)-Gamma(3,2);
a = zeros(4,1);
a(1,1) = .0115;
a(2,1) = .0183;
a(3,1) = .0098;
a(4,1) = .0146;
R = zeros(4,4);
R(1,1) = .8472;
R(1,2) = .1443;
R(1,3) = -.0298;
R(1,4) = -.0669;
R(2,1) = .0431;
R(2,2) = .8774;
R(2,3) = -.0581;
R(2,4) = -.0552;
R(3,1) = .1459;
R(3,2) = .2044;
R(3,3) = .7217;
R(3,4) = -.1515;
R(4,1) = .0177;
R(4,2) = .1861;
R(4,3) = -.0786;
R(4,4) = .7811;
Sigmav = zeros(4,4);
Sigmav(1,1) = .0132;
Sigmav(2,2) = .0154;
Sigmav(3,3) = .0182;
Sigmav(4,4) = .0127;
H = zeros(3,3);
H(1,1) = .9;
H(1,2) = .2;
H(1,3) = .4;
H(2,1) = H(1,2);
H(2,2) = .6;
H(2,3) = .1;
H(3,1) = H(1,3);
H(3,2) = H(2,3);
H(3,3) = .7;
G = zeros(3,3);
G(1,1) = 4;
G(1,2) = .3;
G(1,3) = .1;
G(2,1) = G(1,2);
G(2,2) = 2;
G(2,3) = .05;
G(3,1) = G(1,3);
G(3,2) = G(2,3);
G(3,3) = .5;

t = 1;
T = 100;  

% Model Solution
% Transform System to Canonical Form 
% x[t+tau] = A * x[t+tau-1] + B * E(x[t+tau+1]|I[t+tau]) 
%            + w[t+tau] 
% w(t+tau] = D0 + D1 * m[t+tau]
% m[t+tau] = a + R * m[t+tau-1] + v[t+tau]
M00 = H+(1+beta)*G;
M10 = G;
M01 = beta*M10;
iota = [1; 1; 1];
theta = 1/(iota'*inv(M00)*iota);
N = eye(3)-theta*iota*iota'*inv(M00);
A = M00\N*M10;
B = M00\N*M01;
D0 = M00\N*H*alpha+theta*inv(M00)*iota;
D1 = M00\N*H*[Gamma-delta*omega0' delta];
% Transform System to Canonical Form
% xcurl[t+tau] = Acurl * xcurl[t+tau-1] 
%                + Bcurl * E(xcurl[t+tau+1]|I[t+tau]) 
%                + wcurl[t+tau] 
% wcurl(t+tau] = D0curl + D1curl * m[t+tau]
% m[t+tau] = a + R * m[t+tau-1] + v[t+tau]
P = [1 0; 0 1; -1 -1];
e = [0; 0; 1]; 
D0curl = (P'*P)\P'*(D0+(A+B-eye(3))*e);
D1curl = (P'*P)\P'*D1;
Acurl = (P'*P)\P'*A*P;
Bcurl = (P'*P)\P'*B*P;

sstar = inv(eye(3)-A-B)*(D0+D1*inv(eye(4)-R)*a);
[dim1,dim2] = size(Acurl);
[dim3,dim4] = size(D1curl);

% Carry Out Recursions
captheta = zeros((T-t+1)*dim1,dim1);
capr = zeros((T-t+1)*dim1,1);
capp = zeros((T-t+1)*dim1,1);
capq = zeros((T-t+1)*dim1,dim4);
caps = zeros((T-t+1)*dim1,1);
captheta(1:dim1,:) = eye(dim1);
capr(1:dim1,:) = D0curl;
capp(1:dim1,:) = zeros(dim1,1);
capq(1:dim1,:) = D1curl;
caps(1:dim1,:) = Bcurl*sstar(1:2,1);
for i = 1:T-t
   captheta(i*dim1+1:(i+1)*dim1,:) = eye(dim1)-Bcurl/ ...
                     captheta((i-1)*dim1+1:i*dim1,:)*Acurl;
   capr(i*dim1+1:(i+1)*dim1,:) = ...
                     captheta(i*dim1+1:(i+1)*dim1,:)\Bcurl ...
                     *capr((i-1)*dim1+1:i*dim1,:) ...
                     +captheta(i*dim1+1:(i+1)*dim1,:)\D0curl;
   capq(i*dim1+1:(i+1)*dim1,:) = ...
                     captheta(i*dim1+1:(i+1)*dim1,:)\Bcurl ...
                     *capq((i-1)*dim1+1:i*dim1,:)*R ...
                     +captheta(i*dim1+1:(i+1)*dim1,:)\D1curl;
   capp(i*dim1+1:(i+1)*dim1,:) = ...
                     captheta(i*dim1+1:(i+1)*dim1,:)\Bcurl ...
                     *capq((i-1)*dim1+1:i*dim1,:)*a ...
                     +captheta(i*dim1+1:(i+1)*dim1,:)\Bcurl ...
                     *capp((i-1)*dim1+1:i*dim1,:);
   caps(i*dim1+1:(i+1)*dim1,:) = ...
                     captheta(i*dim1+1:(i+1)*dim1,:)\Bcurl ...
                     *caps((i-1)*dim1+1:i*dim1,:);
end

% Model Simulation
randn('state',20742)
v = randn(4,T-t+1+100);
v(1,:) = Sigmav(1,1)*v(1,:);
v(2,:) = Sigmav(2,2)*v(2,:);
v(3,:) = Sigmav(3,3)*v(3,:);
v(4,:) = Sigmav(4,4)*v(4,:);
m = zeros(4,T-t+1+100);
m(:,1) = inv(eye(4)-R)*a+v(:,1);
for i = 2:T-t+1+100
   m(:,i) = a+R*m(:,i-1)+v(:,i);
end
m = m(:,101:T-t+1+100);
as = zeros(3,T-t+1);
as(1:2,1) = captheta((T-t)*dim1+1:(T-t+1)*dim1,:) ...
            \Acurl*sstar(1:2,1) ...
            +caps((T-t)*dim1+1:(T-t+1)*dim1,:) ...
            +capr((T-t)*dim1+1:(T-t+1)*dim1,:) ...
            +capp((T-t)*dim1+1:(T-t+1)*dim1,:) ...
            +capq((T-t)*dim1+1:(T-t+1)*dim1,:)*m(:,1);          
as(3,1) = 1-as(1,1)-as(2,1);
for i = 2:T-t+1
   as(1:2,i) = captheta((T-t+1-i)*dim1+1:(T-t+2-i)*dim1,:) ...
               \Acurl*as(1:2,i-1) ...
               +caps((T-t+1-i)*dim1+1:(T-t+2-i)*dim1,:) ...
               +capr((T-t+1-i)*dim1+1:(T-t+2-i)*dim1,:) ...
               +capp((T-t+1-i)*dim1+1:(T-t+2-i)*dim1,:) ...
               +capq((T-t+1-i)*dim1+1:(T-t+2-i)*dim1,:)*m(:,i);
   as(3,i) = 1-as(1,i)-as(2,i);          
end
oes = as';

% Display Results
disp('   ')
disp('The optimal expenditure shares are given by: ')
disp('(The i-th expenditure share is in the i-th column, i = 1,2,3, ')
disp('with each column running from t to T.) '), oes
