% =========================================================================
% Simulation of a VAR(p) process
% Impulse Responses and Variance Decomposition analysis

% Fernando Pérez Forero
% fernandojose.perez@upf.edu

% Date: April 19th, 2012

% VAR(p) model
% Y(t)=v+A1*Y(t-1)+...+Ap*Y(t-p)+U(t)
%
% U(t) is N(0,Sigma_u)

% These codes may be freely reproduced only for educational and
% non-comercial research purposes.
% You can use the code or parts thereof for your scientific work as long as
% you acknowledge the corresponding author.

% =========================================================================
clear all
clc

% set the random seed
seed=1950; rand( 'twister' , seed ); randn('state',seed ); 

%% 1. Data Simulation
% ====================

p=4;    % Lag-order
T=100;  % Observations
K=4;    % Variables

% Var-cov matrix
% --------------
scale=2*rand(K,1)+1e-30;
Sigma_U=0.5*rand(K,K)+diag(scale);
Sigma_U=Sigma_U'*Sigma_U;

if all(eig((Sigma_U+Sigma_U')/2)) >= 0
    disp('Sigma_U is positive definite')
else
    error('Sigma_U must be positive definite')
end

Alag=zeros(K,K,p);

lambda=2.5;

for k=1:p
    Alag(:,:,k)=(lambda^(-k))*rand(K,K)-((2*lambda)^(-k))*ones(K,K);
end

Ylag=zeros(K,p);
for k=1:p
    Ylag(:,k)=mvnrnd(zeros(K,1),Sigma_U)';
end

Y=Ylag;
v=ones(K,1);

% Data
for i=1:T    
    temp=v+mvnrnd(zeros(K,1),Sigma_U)';
    for k=1:p
        temp=temp+Alag(:,:,k)*Ylag(:,k);
    end
    
    Y=[Y,temp];
    
    Ylag(:,p)=temp;    
    
    for k=1:p-1
        Ylag(:,k)=Ylag(:,k+1);
    end    
end
clear Ylag temp

A_coeff=v;

for k=1:p
    A_coeff=[A_coeff,Alag(:,:,k)];
end

Y=Y(:,p+1:end);

FS=15;
LW=2;

gr_size2=ceil(K/3);

figure(1)
set(0,'DefaultAxesColorOrder',[0 0 1],...
        'DefaultAxesLineStyleOrder','-|-|-')
set(gcf,'Color',[1 1 1])
set(gcf,'defaultaxesfontsize',FS)
for k=1:K
    subplot(gr_size2,2,k)
    plot(Y(k,:),'b','Linewidth',LW)
    title(sprintf('Y_%d',k))
end

% Companion form
% --------------
F1=[A_coeff(:,2:end)];
Q1=Sigma_U;

if p>1
    % Matrix F - Hamilton eq .[10.1.10]
    F2=cell(1,p-1);
    [F2{:}]=deal(sparse(eye(K)));
    F2=blkdiag(F2{:});
    F2=full(F2);
    F3=zeros(K*(p-1),K);
    F=[F1;F2,F3];    
    
    % Var_cov Q - Hamilton eq .[10.1.11]
    Q2=cell(1,p-1);
    [Q2{:}]=deal(sparse(zeros(K,K)));
    Q2=blkdiag(Q2{:});
    Q2=full(Q2);    
    Q=blkdiag(Q1,Q2);    
    clear F1 F2 F3
else
    F=F1;
    Q=Q1;
    clear F1 Q1
end

if max(abs(eig(F)))<1   % Hamilton eq. [10.1.13]
    disp('VAR is stationary')
else
    error('VAR is not stationary!')
end

% Companion form's Variance-Covariance matrix - Hamilton eq. [10.2.13]

Cap_Sigma0=eye(size(F));
Cap_Sigma1=F*Cap_Sigma0*F'+Q;
diff=max(max(abs(Cap_Sigma1-Cap_Sigma0)));
diff1=[];
epsilon=1e-10;
iter=0;
figure(1)
set(gcf,'Color',[1 1 1])
tic
disp('Iterating covariance matrix ...')

while diff>epsilon    
    iter=iter+1;
    Cap_Sigma0=Cap_Sigma1;
    Cap_Sigma1=F*Cap_Sigma0*F'+Q;
    diff=max(max(abs(Cap_Sigma1-Cap_Sigma0)));
    diff1=[diff1,diff];
    figure(2)
    plot(diff1)
    title('Evolution of convergence','FontSize',14)
end

sprintf('Convergence achieved after %d iterations',iter)

% Alternatively, solve the discrete Lyapunov equation

Cap_Sigma1l=dlyap(F,Q);

% Compare results
diff2=max(max(abs(Cap_Sigma1l-Cap_Sigma1)));

toc

%% 2. Impulse response functions
%  =============================
h=12; % Horizon
P=chol(Sigma_U); % Cholesky decomposition
P=P'; % Following the notation of Lutkepohl - Section 3.7 : U=P*E

capJ=zeros(K,K*p);
capJ(1:K,1:K)=eye(K);

Iresp=zeros(K,K,h);
Iresp(:,:,1)=P;

for j=1:h-1    
    temp=(F^j);
    Iresp(:,:,j+1)=capJ*temp*capJ'*P;    
end
clear temp

FS=15;
LW=2;

% Plot IRF
% --------
figure(3)
set(0,'DefaultAxesColorOrder',[0 0 1],...
        'DefaultAxesLineStyleOrder','-|-|-')
set(gcf,'Color',[1 1 1])
set(gcf,'defaultaxesfontsize',FS-5)
for i=1:K
    for j=1:K
       subplot(K,K,(j-1)*K+i)
       plot(squeeze(Iresp(i,j,:)),'Linewidth',LW)
       hold on
       plot(zeros(1,h),'k')
       title(sprintf('Y_{%d} to U_{%d}',i,j))
       xlim([1 h])
       set(gca,'XTick',0:ceil(h/4):h)
    end
end
ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    text(0.5, 1,'\bf Impulse responses','HorizontalAlignment','center','VerticalAlignment', 'top','FontSize',FS)


%% 3. Variance decomposition
% ==========================

MSE=zeros(K,h);
CONTR=zeros(K,K,h);

for j=1:h   
    % Compute MSE
    temp2=eye(K);
    for i=1:K
        if j==1
            MSE(i,j)=temp2(:,i)'*Iresp(:,:,j)*Sigma_U*Iresp(:,:,j)'*temp2(:,i);
            for k=1:K
            CONTR(i,k,j)=(temp2(:,i)'*Iresp(:,:,j)*P*temp2(:,k))^2;
            end
        else
            MSE(i,j)=MSE(i,j-1)+temp2(:,i)'*Iresp(:,:,j)*Sigma_U*Iresp(:,:,j)'*temp2(:,i);
            for k=1:K
            CONTR(i,k,j)=CONTR(i,k,j-1)+(temp2(:,i)'*Iresp(:,:,j)*P*temp2(:,k))^2; % Lutkepohl eq. (2.3.36)
            end
        end
    end
end

VD=zeros(K,h,K);

for k=1:K
    VD(:,:,k)=squeeze(CONTR(:,k,:))./MSE;  % Lutkepohl eq. (2.3.37)
end

clear temp2 CONTR MSE

% Plot VD
% -------
figure(4)
set(0,'DefaultAxesColorOrder',[0 0 1],...
        'DefaultAxesLineStyleOrder','-|-|-')
set(gcf,'Color',[1 1 1])
set(gcf,'defaultaxesfontsize',FS-5)
for i=1:K
    for j=1:K
       subplot(K,K,(j-1)*K+i)
       plot(squeeze(VD(i,:,j)),'Linewidth',LW)
       hold on
       plot(zeros(1,h),'k')
       title(sprintf('U{%d} to Var(e_{%d,t+h}).',j,i))
       xlim([1 h])
       set(gca,'XTick',0:ceil(h/4):h)
       ylim([0 1])
       set(gca,'YTick',0:0.25:1)
    end
end
ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off');
    text(0.5, 1,'\bf Variance Decomposition','HorizontalAlignment','center','VerticalAlignment', 'top','FontSize',FS)