function [MOM]=model(THETA,NSTR);

global nobs init slong nsim nsec poids trunc nvar Sel XT DES CRE NTM 
global eu ew;
global mom_h Omega
global r s kom phi0 nss beta wb rhow sew;

s_str=size(NSTR,1);
for i=1:s_str;
   eval([NSTR(i,:),'= THETA(',int2str(i),');']);
end;

abar	= 0;%sqrt(abar*abar);
rhou	= sqrt(rhou*rhou);
alp1	= sqrt(alp1*alp1);
phi1	= sqrt(phi1*phi1);
s		= sqrt(s*s);
seu	= sqrt(seu*seu);
kom	= sqrt(kom*kom);
a2		= abar/(r+s);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                              Modèle Théorique                             %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I		= [0:trunc]';
I1		= I+1;
T		= 1:slong;
b0		= (beta*(1-s)).^I;
bs		= b0.*exp(0.5*seu*seu*(1-rhou.^(2*I1))/(1-rhou*rhou));
bw		= b0.*exp(0.5*sew*sew*(1-rhow.^(2*I1))/(1-rhow*rhow));
er		= rhou.^I1;
es		= rhow.^I1;

CT		= zeros(slong+1,1);   
DT		= zeros(slong+1,1);
NT		= zeros(slong+1,1);
N		= zeros(slong+1,nsec);
XT		= zeros(slong+1,nsec);
ui		= 0;
uw		= 0;
N(1,:)= (nss/nsec)*ones(1,nsec);
for i	= 1:slong+1;
   ui			= rhou*ui+seu*eu(i,:);
   uw			= rhow*uw+sew*ew(i);
   %xt = a0*ui+a2;
   xt			= beta*((bs'*exp(er*ui))-wb*(bw'*exp(es*uw)))+a2;
   XT(i,:)	= xt;
   NT(i)		= N(i,:)*poids;
   NE			= (xt>0).*(kom*(xt.^(alp1/(1-alp1)))*(1-NT(i)));
   NL			= -(xt<-phi0).*((1-s)*(1-(1-(xt+phi0)/phi1).^(-0.5))).*N(i,:);
%	NL			= -(xt<-phi0).*((1-s)*N(i,:).*((-xt).^(1/phi1))./(1+((-xt).^(1/phi1))));
   N(i+1,:)	= (1-s).*N(i,:)+NE+NL;
   CT(i)		= ((xt>0).*(kom*(xt.^(alp1/(1-alp1))))*(1-NT(i)))*poids;
   DT(i)		= ((1-s)*(1-(1-(xt+phi0)/phi1).^(-0.5)).*(xt<-phi0).*N(i,:)+s*N(i,:))*poids;
%   DT(i)		= ((1-s)*N(i,:).*((-xt).^(1/phi1))./(1+((-xt).^(1/phi1))).*(xt<-phi0)+s*N(i,:))*poids;
end;
%
% Calcul des variables agrégées
%
NT0	= NT(init:slong-1,:); 
NT1	= NT(init+1:slong,:);
NTM	= 0.5*(NT0(:)+NT1(:));
%NTM	= NT0(:);
CRE	= 100*CT(init:slong-1)./NTM;
DES	= 100*DT(init:slong-1)./NTM;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                           %
%                              Moments simulés                              %
%                                                                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lt=length(CRE);

%1-2
m1cre=mean(CRE);
m1des=mean(DES);

%CRE=CRE-m1cre;
%DES=DES-m1des;

C2=CRE.*CRE;
C3=C2.*CRE;
C4=C3.*CRE;
D2=DES.*DES;
D3=D2.*DES;
D4=D3.*DES;

%3-4
m2cre=mean(C2);
m2des=mean(D2);
%5-6
m3cre=mean(C3);
m3des=mean(D3);
%7-8
m4cre=mean(C4);
m4des=mean(D4);
%9-10
r1cre=mean(CRE(2:lt).*CRE(1:lt-1));
r1des=mean(DES(2:lt).*DES(1:lt-1));
%11-12
r2cre=mean(C2(2:lt).*C2(1:lt-1));
r2des=mean(D2(2:lt).*D2(1:lt-1));
%13-15
c1pn=mean(CRE.*DES);
c1pnr=mean(CRE(2:lt).*DES(1:lt-1));
c1pna=mean(CRE(1:lt-1).*DES(2:lt));
%16-18
c2pn=mean(C2.*DES);
c2pnr=mean(C2(2:lt).*DES(1:lt-1));
c2pna=mean(C2(1:lt-1).*DES(2:lt));
%19-21
c3pn=mean(D2.*CRE);
c3pnr=mean(D2(2:lt).*CRE(1:lt-1));
c3pna=mean(D2(1:lt-1).*CRE(2:lt));
%22-24
c4pn=mean(C2.*D2);
c4pnr=mean(C2(2:lt).*D2(1:lt-1));
c4pna=mean(C2(1:lt-1).*D2(2:lt));

MOM=[m1cre;m1des;
   m2cre;m2des;
   m3cre;m3des;
   m4cre;m4des;
   r1cre;r1des;
   r2cre;r2des;
   c1pn;c1pnr;c1pna;
   c2pn;c2pnr;c2pna;
   c3pn;c3pnr;c3pna;
   c4pn;c4pnr;c4pna];

MOM=MOM(Sel);