
function [irf,irf1,irf2,irf_all]=run_LP_MR(data0,id,T0,id_m,ord,shock,normvar,names,varargin)


y=data0(:,id);
REPS=10000;
BURN=1000;
SKIP=1;
%-MODEL SPECIFICATION-----------------------------------------------------%

%declare shock variable & shock size
shockSize =1;           %percentage points
N=length(id);
%set VAR specification, IRFs horizon, 
modelSpec.nVARlags      =12; %number of lags in VAR (& prior)
modelSpec.nBLPlags      =12; %number of lags in BLP
modelSpec.nHorizons     =36;
modelSpec.bandsCoverage =90;
%set BLP prior type
modelSpec.priorType ='VAR'; %BLP only, alternatives: 'VAR', 'RW'; 
%initialize prior on presample
modelSpec.presample=true;

%-IDENTIFICATION SCHEME---------------------------------------------------%
%choose identification scheme
modelSpec.identification      ='PSVAR'; %alternatives: 'PSVAR', 'CHOL';
if~isempty(varargin)
    proxytype=varargin{1};
    modelSpec.proxytype=proxytype;
end
%if PSVAR choose external instrument
modelSpec.m_t  =data0(T0+1:end,id_m); 
modelSpec.ord=ord;
modelSpec.shock=shock;
%-ASSEMBLE DATA-----------------------------------------------------------%
dataStructure.data=y(T0+1:end,:);
dataStructure.varname=names(id);
dataStructure.preSdata=y(1:T0,:);
tmp=zeros(1,N);
tmp(shock)=1;
selmat=logical(tmp);


%load shock variable and normalization in model structure 
modelSpec.shockSize     =shockSize*selmat';
modelSpec.shockVar      =selmat;
modelSpec.normvar=normvar;
%-HYPERPRIORS SETTINGS----------------------------------------------------%

%random walk prior
hyperPars.isrw      =true(1,length(dataStructure.varname)); %rw prior

%hyperpriors initial values
%
hyperPars.lambda    =.1;    %tightness of VAR coeffs prior (the higher lambda the closer to OLS)
hyperPars.lambdaC   =1e5;   %intercept (you want this to be large)
hyperPars.lambdaP   =.4;    %tightness of projCoeffs prior (same role as above)
%
hyperPars.miu       =[];     %sum of coefficients prior (constraint multiplier)
hyperPars.theta     =[];     %cointegration prior (constraint multiplier) 
hyperPars.alpha     =[];     %lag decaying coeff for NIW prior


%set hyperpriors options (matches GLP fields); if you want default values
%(when available) set to empty [];
hyperPriorsOptions.hyperpriors   = true;                  %find optimal hyperparameters: NO default option
hyperPriorsOptions.Vc            = 1e5;                   %variance of the VAR constant (default=1e6)
hyperPriorsOptions.pos           = find(~hyperPars.isrw); %position of stationary variables
hyperPriorsOptions.MNalpha       = [];                    %lag decaying coeff of NIW prior (default=2)
hyperPriorsOptions.MNpsi         = false;                 %residual variance univariate AR(1) std (default=hyperprior)
hyperPriorsOptions.noc           = false;                 %sum of coefficients prior: NO default option
hyperPriorsOptions.sur           = false;                 %cointegration prior: NO default option
hyperPriorsOptions.Fcast         = false;                 %build forecasts: NO default option
hyperPriorsOptions.hz            = modelSpec.nHorizons;   %max forecast horizon: NO default option
hyperPriorsOptions.mcmc          = false;                 %run metropolis-hasting algorithm: NO default option
hyperPriorsOptions.Ndraws        = REPS;                  %default=20k
hyperPriorsOptions.Ndrawsdiscard = BURN;                   %default=10k
hyperPriorsOptions.MCMCconst     = SKIP;                     %default=1
hyperPriorsOptions.MCMCfcast     = false;                 %store forecast at each MCMC draw (default=true)
hyperPriorsOptions.MCMCstorecoeff= false;                 %store coefficients at each MCMC draw (default=true)
hyperPriorsOptions.initialValues = hyperPars;             %see above


%sampling from paramenters distribution
GibbsOptions.iterations          = REPS;
GibbsOptions.burnin              = BURN;
GibbsOptions.jump                = SKIP;

hyperPriorsOptions.GibbsOptions  = GibbsOptions;



%-------------------------------------------------------------------------%
%-------------------------------------------------------------------------%
%-------------------------------------------------------------------------%

modelSpec.dataStructure=dataStructure;

dataMatrix = dataStructure.data; 


%-ESTIMATE IRFs-----------------------------------------------------------%


%1) BLP with conjugate prior (BLP)

irfBLP=IRFbayesianLocalProj_5(dataMatrix,modelSpec,hyperPriorsOptions);

irf=irfBLP.irfs;
irf1=irfBLP.irfs_l;
irf2=irfBLP.irfs_u;

irf_all(1,:,:)=irf(:,end);
irf_all(2,:,:)=irf1(:,end);
irf_all(3,:,:)=irf2(:,end);

