%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                    Déclaration des variables globales                     %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global nobs init slong nsim nsec poids trunc nvar XT DES CRE NTM 
global eu ew;
global mom_h Omega
global r s kom phi0 nss beta wb rhow sew Sel abar;
ficm	= 'model';			% Nom du fichier de modèle
fic	= 'sgmm.out';		% Nom du fichier de sortie
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                    Chargement des moments sur les données                 %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%usflows3;
usflows2;
Sel	= [1:24];
%Sel	= [1:8 13];
mom_h	= mom_h(Sel);
%Omega	= inv(diag(diag(mv_h(Sel,Sel)))); 
Omega	= inv(mv_h(Sel,Sel)); 
dim_j	= size(mom_h,1);

%Omega=diag(diag(ones(4)));
%
% Nombre de simulations
%
nobs	= 87;					% # of obs.
init	= 250;				% # initial conditions
nsim	= 10;				% # of sim.
slong	= init+nsim*nobs;	% # total length
nsec	= 2;					% # of sectors
trunc	= 1000;				% # truncation of the series
poids	= ones(nsec,1)/nsec;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                    Déclaration des paramètres initiaux                    %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nss	= 0.57;
r		= 0.01;
beta	= 1/(1+r);
wb		= 1;
phi0	= 0;
rhow	= 0.852398;
sew	= 0.0072878;

rhou	= 0.7965;
seu	= 0.1966;
s		= 0.0480;
alp1	= 0.2498;
phi1	= 8.4453;
abar	= 0;
kom	= 0.3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                    	   Tirage des innovations                           %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
randn('state',2);
e	= randn(slong+2,2);
for i=1:2;
   e(:,i)=(e(:,i)-mean(e(:,i)))/std(e(:,i));
end;
ew	= e(:,1);
eu	= [e(:,2) -e(:,2)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                    Déclaration des paramètres initiaux                    %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

theta0=[
0.29577516850951   
0.15588923271566   
0.05218945737242   
0.41944681241130
1.63597463052307   
];

rhou	= 0.86416114;
seu	= 0.10674130;
alp1	= 0.34555853;
phi1	= 11.91356611;

rhou	= 0.5146;
seu	= 0.1154;
s		= 0.0490;
alp1	= 0.3340;
phi1	= 1.0952;

rhou	= 0.765;
seu	= 0.016;
s		= 0.030;
alp1	= 0.2498;
phi1	= 8.4453;

theta0= [rhou;seu;s;alp1;phi1];
theta0=[
   0.8358508652211
   0.02453210336337
   0.04926111758973
   0.39105337252842
   1.50611924138084
];
NSTRU		= [   
   'rhou  ';
   'seu   ';
   's     ';
   'alp1  ';
   'phi1  ';
%   'abar  ';
%   'kom   ';
];
dim_coef	= size(theta0,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                     Appel de la procédure de minimisation                 %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
opt=foptions;
opt(3)	= 1e-2;
opt(4)	= 1e-4;
opt(14)	= 1e9;
%sgmmj(theta0,2,'model',NSTRU)
%theta=theta0;
theta		= fmins('sgmmj',theta0,opt,[],1,ficm,NSTRU);
%theta		= fminu('sgmmj',theta0,opt,[],1,ficm,NSTRU);
nn_par	= dim_coef;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                          Calcul des Ecarts-type                           %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%dJ			= numgrad('sgmmj',theta,'c',1,ficm,NSTRU);

mom_m		= sgmmj(theta,2,ficm,NSTRU);
df			= numgrad('sgmmj',theta,'c',2,ficm,NSTRU);
dig		= df';
%
% Statistiques pour les paramètres structurels
%
mvc_str	= inv(dig*Omega*dig');
std_str	= sqrt(real(diag(mvc_str)))';
T_str		= real(theta(:)./std_str(:));
%
% Statistiques pour les paramètres auxiliaires engendrés par le modèle
%
mvc_m		= real(dig'*mvc_str*dig);
std_m		= sqrt(diag(mvc_m));
T_m		= real(mom_m(:)./std_m(:));
%
% Statistiques pour les paramètres auxiliaires sur les données
%
std_h		= sqrt(diag(inv(Omega)));
T_h		= real(mom_h(:)./std_h(:));
%
% Test de diagnostique
%
VD_hm		= inv(Omega)-mvc_m;
std_diag	= sqrt(diag(VD_hm));
diff_mom	= (mom_m(:)-mom_h(:));
T_diag	= diff_mom(:)./std_diag(:);
%
% Test de sur-identification globale
%
J_stat	= (nsim/(1+nsim))*real(diff_mom'*Omega*diff_mom); 
disp(''); 
disp('');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                          Affichage des résultats                          %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%eval(['delete ',fic])
eval(['diary ',fic]);
clc
fprintf('\n\t\tPARAMETRES STRUCTURELS \n\n');
disp(repmat('-',1,60))
fprintf('Coeff.\t\tValeur\t\tEcart-Type\tT-Stat.\n')
disp(repmat('-',1,60))
for i=1:dim_coef;
   disp(sprintf('%s\t:\t%8.4f\t%8.4f\t%8.4f',NSTRU(i,:),theta(i),std_str(i),T_str(i)))
end
disp(repmat('-',1,60))
fprintf('J-stat : %10.4f\t degré de lib. : %d\t\n',J_stat,dim_j-dim_coef)
disp(repmat('-',1,60))
fprintf('\n\n\t\tPARAMETRES AUXILIAIRES\n\n');
fprintf('\n\n\t\t  DANS LES DONNEES\n\n');
fprintf('Valeur\t\tEcart-Type\tT-Stat.\n')
disp(repmat('-',1,60))
for i=1:dim_j;
   disp(sprintf('%8.4f\t%8.4f\t%8.4f',mom_h(i),std_h(i),T_h(i)))
end
disp(repmat('-',1,60))
fprintf('\n\n\t\t  DANS LE MODELE\n\n');
fprintf('Valeur\t\tEcart-Type\tT-Stat.\n')
disp(repmat('-',1,60))
for i=1:dim_j;
   disp(sprintf('%8.4f\t%8.4f\t%8.4f',mom_m(i),std_m(i),T_m(i)))
end
disp(repmat('-',1,60))

disp('Test de diagnostique')
disp(repmat('-',1,60))
for i=1:dim_j;
   disp(sprintf('%8.4f',T_diag(i)))
end
disp(repmat('-',1,60))
disp(sprintf('Rang du jacobien : %d',rank(dig)))
diary off

diary off

disp(sprintf('rhou\t= %6.8f;',theta(1)));
disp(sprintf('seu\t= %6.8f;',theta(2)));
disp(sprintf('s\t= %6.8f;',theta(3)));
disp(sprintf('alp1\t= %6.8f;',theta(4)));
disp(sprintf('phi1\t= %6.8f;',theta(5)));
disp(sprintf('abar\t= %6.8f;',theta(6)));
%disp(sprintf('kom\t= %6.8f;',theta(7)));
