% ARRW_IMPORTANCE_SAMPLER0.M
% Implements importance sampler for models that satisfy zero restrictions
% and sign restrictions as in Arias et al. (2018)

clear;

% Load all draws that satisfy these restrictions
load posterior0

% ARRW notation
nlag      = p;              % number of lags
nvar      = n;               % number of endogenous variables
m         = nvar*nlag + 1;   % number of exogenous variables

% Define Z(1) like E(1) on p. 480 of Kilian and Luetkepohl (2017)
Z=cell(nvar,1);
for i=1:nvar
    Z{i}=zeros(0,numel(0)*nvar);  
end;
nz1=2; % Number of zero restriction for shock 1
Z{1}=zeros(nz1,numel(0)*nvar);
Z{1}(1,1)=1;
Z{1}(2,2)=1;

% Number of zero restrictions
nzeros=0;
for i=1:nvar
    nzeros=nzeros+size(Z{i},1);
end;

% computes the dimension of space in which the spheres live that map into
% the orthogonal matrices satisfying the restriction.  
dim=0;
for i=1:nvar
    dim=dim+nvar-(i-1+size(Z{i},1));
end

% Get random W 
W=cell(nvar,1);
for j=1:nvar
    W{j}=randn(nvar-(j-1+size(Z{j},1)),nvar);
end;

% Function handles
fo=@(x)f_h(x,nvar,m);
fo_str2irfs_inv=@(x)IRFToStructural(x,nvar,nlag);
fs=@(x)ff_h(x,nvar,m,nlag,W,dim,Z); 
r=@(x)ZeroRestrictions(x,nvar,nzeros,m,nlag,Z)';
r_irfs=@(x)IRFRestrictions(x,Z); 

for j=1:size(z,1)
     
    disp(j)
    
    % IK notation for jth model
    A      = (reshape(z(j,1:(1+n*p)*n),1+n*p,n))'; 
    B0inv  = reshape(z(j,(1+n*p)*n+n*n+1:(1+n*p)*n+2*n*n),n,n);    
    A=[A(:,2:end) A(:,1)]; % RF level parameters rearranged with constant ordered last
    
    % ARRW notation:
    % yt(t) A0 = xt(t) Aplus + constant + et(t) for t=1...,T;
    % yt(t)    = xt(t) B     + ut(t)            for t=1...,T;
    % x(t)     = [yt(t-1), ... , yt(t-nlag), constant];
    
    % Translating IK notation to ARRW notation for structural parameters A0 and Aplus 
    A0=inv(B0inv)';   % Tranpose required since ARRW use row format
    Aplus=[];
    for i=1:nlag
        Aplus=[Aplus; (inv(B0inv)*A(:,(i-1)*nvar+1:i*nvar))'];
    end;   
    Aplus=[Aplus; (inv(B0inv)*A(:,end))'];
   
    % Inputs for importance sampler
    structpara=[vec(A0); vec(Aplus)];
    irfpara=StructuralToIRF(structpara,nvar,nlag);

    % Computing importance weights
    % The function LogVolumeElement has inputs: (function,data,function)
    storevephi(j,1)   = LogVolumeElement(fo,structpara)   + LogVolumeElement(fo_str2irfs_inv,irfpara);        
    storevegphiZ(j,1) = LogVolumeElement(fs,structpara,r) + LogVolumeElement(fo_str2irfs_inv,irfpara,r_irfs); 
    uw(j,1)           = exp(storevephi(j,1) - storevegphiZ(j,1));
      
end;

% Compute the normalized weights 
imp_w  = uw/sum(uw);

% Importance resampler
za=[]; 
for s=1:size(z,1)
    is_draw  = randsample(1:length(imp_w),1,true,imp_w);  % is_draw is scalar index number 
    za=[za; z(is_draw,:)];
end;

% Save resampled results
save posterior_arrw0 M rep p n H t Ydep X za; 
 
