%% VECM_demo.mlx
% Joel Hasbrouck
% November 7, 2018
% Simulate a small number of cointegrated prices
% Estimate a VECM of order P using Matlab's Econometrics Toolbox vecm routines.
% Construct a random-walk decomposition from the Matlab vecm estimates.
% Estimate the same VECM using the sparse routines in an MVARi object.
% Construct a random-walk decomposition from the MVARi estimates.
%% Initializations and simulations

clc; clear all; close all
format compact
addpath('./mClasses','./mFiles')
rng('default')
mv = MVARi;
T=5000;
spd.setgetMax(T);
nPrices = 3;
density = .6;
mv.sim(density,nPrices,0,5) % The first price is the leader; other obs are randomly lagged between 0 and 3 intervals.
mv.pricePlot

%
%% Estimation
% Estimate with matlab VECM routines

P1 = 3; % order of AR
first = 100-P1; last=4900;  % defines window over which data will be analyzed.
xs = [];
for i=1:mv.nPrices
    xs=[xs mv.prices(i).toCol];
end
xs = xs(first:last,:);
xs = xs - repmat(mean(xs),size(xs,1),1);    % demean the data.
n = size(xs,2);
%   set up model
mdl = vecm(n,mv.nPrices-1,P1);
B = [ones(n-1,1) -eye(n-1)];
mdl.Cointegration = B';
mdl.Constant = zeros(n,1);
mdl.Trend = zeros(n,1);
mdl.SeriesNames = mv.priceNames;
[est,se] = estimate(mdl,xs,'Display','full','Model','H2','MaxIterations',1);
% [est,se] = estimate(mdl,xs,'Display','full','MaxIterations',1);

% Now perform the same estimation using the sparse routines.
% The sparse routines are set up using polynomial distributed lag specifications. Here, we want the usual AR terms. This is done by setting the polynomials to that the the design matrix is the identity matrix of order P1.

pdls = polynom.empty();
for k=1:P1
    pdls = horzcat(pdls,polynom(0,1,['pdl' int2str(k)],k-1));
end

% The design matrix works out to the identity matrix:

pdls.designMatrix

% Now finish the estimates.

mv.polys = pdls.copy();
mv.intercept = false;
mv.ecm = true;
mv.setup(true);
mv.setNamesPDL(true);
mv.zpzLayoutPDL();
mv.firstValid = first; mv.lastValid = last;
mv.eVecDemean(true);
zpzDisplay = 1;
mv.buildzpzSparsePDLpar(zpzDisplay,0);
mv.estimateSparse();
mv.dispEstimates

% The two sets of estimates are quite close, but not identical to numerical precision. I have not determined the source of the discrepancies.
%
%% Random-walk decompositions
% The random-walk decomposition for the Toolbox vecm:

r = randomWalkDecomp;
r.initm(est);
r.isBoundsGrouped({{'p*'},{'p2','p3'}});
r.rDisplay('rwd from the toolbox vecm');

% The random-walk decomposition for the sparse (MVARi) vecm

r=randomWalkDecomp;
r.init(mv);
r.isBoundsGrouped({{'p*'},{'p2','p3'}});
r.rDisplay('rwd from the sparse vecm');

%
%% IRF Computation (MVARi only)

nAhead = 20;
mv.irfPDLpacked(nAhead, 1, 0);
irf=squeeze(mv.irfPacked(3,:,1));
irft=mv.irft;
plot(irft,irf)
title('p3 after a one-unit shock to p*')
