var v1 v2;
varexo eps1 eps2 eps3;

parameters rho11 rho12 sigma1 sigma2 phi rho;

% if    
%   v1 = rho11*v1(-1) + rho12*v2(-1) + sigma1*eps1 + phi*eps3;
%   v2 = rho12*v1(-1) + rho11*v2(-1) + sigma2*eps2 - phi*eps3;
% standard deviation of shocks 
%   sqrt( sigmai^2+phi^2 )
% covariance of the two shocks 
%      E[ (sigma1*eps1 + phi*eps3) * (sigma2*eps2 - phi*eps3) ] = -phi^2
% correlation of the two shocks
%     -phi^2 / ( sqrt(sigma1^2+phi^2) * sqrt(sigma2^2+phi^2) )
%   = -phi^2 / ( sigma1*sigma2+phi^2 )
% since
%  rho  = -phi^2 / (sigma1*sigma2+phi^2)
% then
%  phi  = sqrt( -sigma1*sigma2*rho/(1+rho) )

rho11  = 0.8;
rho12  = 0.1;
sigma1 = 0.3;
sigma2 = 0.3;
phi    = 0.2; 
rho    =-phi^2 / (sigma1*sigma2+phi^2);

stdv_shocks = sqrt(sigma1*sigma2+phi^2);
corr_shocks =-phi^2/(sigma1*sigma2+phi^2);

model;
   #phi1 = sqrt( -sigma1*sigma2*rho/(1+rho) );
   v1 = rho11*v1(-1) + rho12*v2(-1) + sigma1*eps1 - phi1*eps3;
   v2 = rho12*v1(-1) + rho11*v2(-1) + sigma2*eps2 - phi1*eps3;
end;

shocks ;
   var eps1 = 1;
   var eps2 = 1;
   var eps3 = 1;
end;

steady;

check;

% simulate and save data
stoch_simul(periods=1000,order=1);
save simdata v1 v2;

% estimation
varobs v1 v2;

% ML estimation
estimated_params ;
   rho11,  0.95;
   rho12,  0.02;
   sigma1, 0.2;
   sigma2, 0.2;
   rho, -0.1;
end; 
%estimation(datafile=simdata);

% Bayesian estimation
estimated_params;
   rho11, beta_pdf, 0.7, 0.1;
   rho12, beta_pdf, 0, 0.05, -1, 1;
   sigma1, inv_gamma_pdf, 0.2, Inf;
   sigma2, inv_gamma_pdf, 0.2, Inf;
   rho, -0.2, beta_pdf, 0, 0.3, -1, 1;
end; 
options_.mh_nblck     = 4;
options_.mh_replic    = 5000;
options_.mh_jscale    = 1.2;
options_.mode_compute = 1;
estimation(datafile=simdata,mode_check);

% identification(advanced=1);
identification(advanced=1,parameter_set = posterior_mean);
% dynare_sensitivity(identification=1, morris=2);