/* ******************************************************************* */
/*   CHAPITRE 19 : DETERMINANTS DES DECISIONS D'OFFRE DE TRAVAIL       */
/* ******************************************************************* */

libname exo 'Z:\Econometrie appliquee\donnees';

ods pdf file="Z:\Econometrie appliquee\Resultats\chapitre19.pdf";
title;

************************************************************************
*                 Intitulé des séries                                  *
************************************************************************

*	l = nombre d'heures de travail réalisées par la femme en 1995      ;
*	enf6 = nombre d'enfants de moins de 6 ans du ménage                ;
*	enf618 = nombre d'enfants entre 6 et 18 ans du ménage              ;
*	agef = âge de la femme                                             ;
*	form = nombre d'années de formation de la femme                    ;
*	w = salaire horaire de le femme (en francs 1995)                   ;
*	revenu = revenu du ménage en 1995                                  ;
*	chomage = taux de chomage dans la région                           ;

* ******************************************************************** ; 
* Problème n°19.1 : Construction des variables et analyse descriptive  ;
* ******************************************************************** ; 

* Point 1 : construction des variables ;
* ------------------------------------ ;
data chapit19;
  set exo.chapit19;
  if L > 0 then LFP=1;
    else LFP=0;
  V = revenu-(l*w);
  if form > 12 then f=1;
    else f=0;
  A2 = agef * agef;
  if LFP=1 then NLFP=0;
    else NLFP=1;
run;

proc freq data=chapit19;
  table lfp f;
run;
 
* Point 2 : calcul des moyennes ;
* ----------------------------- ;
proc means data=chapit19;
run;

proc means data=chapit19;
  class lfp;
run;


* ******************************************************************** ; 
* Problème n°19.2 : Modèle de participation                            ;
* ******************************************************************** ; 

* Point 1 : Modèle dichotomique logit et probit ;
* --------------------------------------------- ;
proc logistic data=chapit19 descending outest=beta1 covout;
  model lfp = enf6 enf618 agef a2 f v chomage / technique=newton
                                                rsquare 
												stb;
  output out=res1 p=pred1;
  title "Modèle de participation avec modélisation de la loi logistique";

run;

proc logistic data=chapit19 descending outest=beta2 covout;
  model lfp = enf6 enf618 agef a2 f v chomage / link=probit
                                                technique=newton
                                                rsquare 
												stb;
  output out=res2 p=pred2;
  title "Modèle de participation avec modélisation de la loi normale";
run;

* Point 2 : Effet global et marginal de l'âge ;
* ------------------------------------------- ;
data beta2b;
  set beta2(firstobs=1 obs=1);
  eage=-agef/a2;
  agemax=-agef/(2*a2);
run;

proc print data=beta2b noobs;
  var agef a2 eage agemax;
  title "Effets global et marginal de l'age sur la probabilité de participer";
run;


* Point 3 : effets marginaux des variables au point moyen ;
* ------------------------------------------------------- ;
proc means data=chapit19 noprint;
  var enf6 enf618 agef a2 v chomage;
  output out=eff_marginal;
run;

data eff_marginal2;
  set eff_marginal(firstobs=4 obs=4);
  drop _TYPE_ _FREQ_;
  rename enf6=m_enf6 enf618=m_enf618 agef=m_agef a2=m_a2 v=m_v chomage=m_chomage;
run;

data eff_marginal3;
  set eff_marginal2; set beta2b;
  drop _LINK_ _TYPE_ _STATUS_ _NAME_ _LNLIKE_ eage agemax;
run;

* Effets marginal de la variable qualitative FORMATION;
* ----------------------------------------------------;
data eff_marginal3;
  set eff_marginal3;
  mxb = intercept + enf6*m_enf6 + enf618*m_enf618 + agef*m_agef + a2*m_a2 + v*m_v + chomage*m_chomage;
  pm=PROBNORM(mxb);
  mxb1 = intercept + enf6*m_enf6 + enf618*m_enf618 + agef*m_agef + a2*m_a2 + v*m_v + chomage*m_chomage + f;
  pm1=PROBNORM(mxb1);
  ff=pm1-pm;
run;

proc print data=eff_marginal3 noobs;
  var mxb mxb1 pm pm1 ff;
  title "Effet marginal de la variable qualitative FORMATION";
run; 

* Effets marginaux des variables quantitatives au point moyen de l'échantillon;
* ----------------------------------------------------------------------------;
data eff_marginal4;
  set eff_marginal3;
  dmxb    = PDF('Normal',mxb);
  emenf6  = enf6*dmxb;
  emenf618= enf618*dmxb;
  emagef  = (agef+a2*m_agef)*dmxb;
  emv     = v*dmxb;
  emc     = chomage*dmxb;
run;

proc print data=eff_marginal4 noobs;
  var emenf6 emenf618 emagef emv emc;
  title "Effets marginaux des variables continues au point moyen de l'échantillon";
run; 


* Quasi-elasticité des variables quantitatives au point moyen de l'échantillon;
* ----------------------------------------------------------------------------;
data eff_marginal4;
  set eff_marginal4;
  qeenf6  = emenf6*m_enf6;
  qeenf618= emenf618*m_enf618;
  qev     = emv*m_v;
  qec     = emc*m_chomage;
run;

proc print data=eff_marginal4 noobs;
  var qeenf6 qeenf618 qev qec;
  title "Quasi-elasticité des variables continues au point moyen de l'échantillon";
run; 

* Effets marginaux des variables quantitatives pour une probabilité donnée de 0.90;
* --------------------------------------------------------------------------------;
data eff_marginal4;
  set eff_marginal4;
  qt     = PROBIT(0.9);
  dqt    = PDF('Normal',qt);
  empenf6  = enf6*dqt;
  empenf618= enf618*dqt;
  empagef  = (agef+a2*m_agef)*dqt;
  empv     = v*dqt;
  empc     = chomage*dqt;
run;

proc print data=eff_marginal4 noobs;
  var empenf6 empenf618 empagef empv empc;
  title "Effets marginaux des variables continues pour une probabilité donnée de 0.90";
run; 

* Quasi-elasticité des variables quantitatives pour une probabilité donnée de 0.90;
* --------------------------------------------------------------------------------;
data eff_marginal4;
  set eff_marginal4;
  qepenf6  = empenf6*m_enf6;
  qepenf618= empenf618*m_enf618;
  qepv     = empv*m_v;
  qepc     = empc*m_chomage;
run;

proc print data=eff_marginal4 noobs;
  var qepenf6 qepenf618 qepv qepc;
  title "Quasi-elasticité des variables continues pour une probabilité donnée de 0.90";
run; 


* Point 4 : Mesure de qualité d'ajustement du modèle ;
* -------------------------------------------------- ;
data res2;
  set res2;
  if pred2 > 0.5 then lfp_predit=1;
    else lfp_predit=0;
 run;

 proc freq data=res2;
   tables lfp_predit*lfp / out=qual nocol norow;
   title "Matrice de confusion";
 run;

data qual2;
  set qual(firstobs=1 obs=1);
  nb1=count;
  fq1=percent;
run;

data qual3;
  set qual(firstobs=4 obs=4);
  nb2=count;
  fq2=percent;
run;

data qual4;
  set qual2; set qual3;
  nb=nb1+nb2;
  fq=fq1+fq2;
  drop count percent lfp lfp_predit;
run;

proc print data=qual4 noobs;
  var nb fq;
  title "Nombre et pourcentage de bonnes prédictions";
run;


* ******************************************************************** ; 
* Problème n°19.3 : Estimation du salaire et de la fonction d'offre    ;
*                   de travail                                         ;
* ******************************************************************** ; 

* Point 1 : estimation du salaire par un modèle Tobit ;
* --------------------------------------------------- ;
proc sort data=chapit19;
  by W;
run;

proc reg data=chapit19;
  model W = agef A2 F / sse mse dw;
  output out=beta4 student=stud p=pred r=e;
  plot student.*predicted.;
  plot r.*W;
  title "Regression par les MCO";
run;
quit;

proc lifereg data=chapit19;
  model W = agef A2 F / d=normal corrb;
  output out=res3 xbeta=beta;
run;

proc qlim data=chapit19 outest=param;
  model W = agef A2 F ;
  endogenous W ~ censored(lb=0);
  output out=result;
run;

data param2;
  set param;
  keep _Sigma sigma Intercept AGEF A2 f;
  sigma=_Sigma*_Sigma;
  rename intercept=b1
         agef=b2
		 a2=b3
		 f=b4;
run;

data param3;
  set param2; set eff_marginal4;
  keep sigma b1 b2 b3 b4 m_agef m_a2;
run;

data param3;
  set param3;
  xsig=(b1+b2*m_agef+b3*m_a2)/sqrt(sigma);
  dsig=probnorm(xsig);
  xsig20=(b1+b2*20+b3*20*20)/sqrt(sigma);
  dsig20=probnorm(xsig20);
  ema42=(b2+2*b3*m_agef)*dsig;
  ema20=(b2+2*b3*20)*dsig20;
run;

proc print data=param3 noobs;
  var ema42 ema20;
  title1 "Effets marginaux pour une femme sans emploi";
  title2 "avec une variable égale à 42 ans puis 20 ans";
run;

* Point 2 : Modélisation selon la méthode d'Heckman ;
* ------------------------------------------------- ;

* Etape 1 : estimation de la propension a accepter le principe ;
* le principe de dedommagement par le modèle PROBIT            ;
* ------------------------------------------------------------ ;
proc sort data=chapit19;
  by descending lfp;
run;

* Estimation du modele par la PROC LOGISTIC ;
* ----------------------------------------- ;
proc logistic data=chapit19 order=data;	
  class lfp;
  model lfp = enf6 enf618 agef a2 f v  / link=probit rsquare technique=newton ;
  output out=imr xbeta=gammaw;
  title "Modélisation Probit de LFP";
run;

proc probit data=chapit19 order=data outest=probit_est;	
  class lfp;
  model lfp = enf6 enf618 agef a2 f v ;
  output out=imr xbeta=gammaw;
  title "Modélisation Probit de LFP";
run;

* Calcul du ratio de Mills ;
* ------------------------ ;
data x(keep=intercep lfp enf6 enf618 agef a2 f v) 	/* variables de l'etape 1 y compris la var_cible */
     w(keep=intercep enf6 enf618 agef a2 f v) 		/* variables de l'etape 1 y compris la var_cible */
     xstar(keep=intercep enf6 enf618 agef a2 v w)  /* variables quanti de l'etape 2 */
     delta(keep=delta)
     h(keep=h)
   /* Dans table B : liste des variables des etapes 1 et 2  */
	 b(keep=L lambda enf6 enf618 agef a2 f v w);
   retain intercep lpf lambda lfp enf6 enf618 agef a2 f v w ;
   set imr;

   /* Calcul du ratio de Mills */
   if (lfp eq 1) then
      lambda=(1/sqrt(2*3.141592654)*exp(-1*gammaw**2/2))/probnorm(gammaw);
     else if (lfp eq 0) then
      lambda=(1/sqrt(2*3.141592654)*exp(-1*gammaw**2/2))/(probnorm(gammaw)-1);
     else lambda=.;

   /* Creation de la constante pour l'utilisation dans les matrices croisées */
   intercep=1;

   /* Creation de h pour l'estimation de la matrice asymétrique de var-covariances */
   /* des coefficients Probit                                                      */
   h=lambda**2+lambda*gammaw;

   /* Creation de delta pour l'estimation de la matrice asymétrique de var-covariances */
   /* des coefficients MCO ; ceci est redondant mais rend les notations suivantes plus */
   /* faciles a utiliser                                                               */ 
   delta=h;

   if (lfp eq 1) then do;	/* Tables SAS selectionnant  */
      output delta;                  	/* les observations pour le  */
      output w;                      	/* calcul des ecart-type du  */
      output xstar;                  	/* modele MCO                */
      output b;
   end;
   output x;                         	/* Tables SAS selectionnant     */
   output h;                         	/* toutes les observationsr     */
                                     	/* pour calculer les ecart-type */
                                     	/* du modele Probit             */
run;

* Etape 2 : estimation du montant par le modèle MCO ;
* ------------------------------------------------- ;
proc reg data=b outest=mco_est;
   model L= lambda enf6 enf618 agef v w  ; 
   output out=err p=pred residual=e;
run;
quit;

proc means data=err maxdec=2; 
  var pred L v w;
  output out=err2;
run;

data err3;
  set err2(firstobs=4 obs=4);
  drop _TYPE_ _FREQ_ _STAT_;
  rename pred=m_pred l=m_l v=m_v w=m_w;
run;

data mco_est;
  set mco_est;
  rename v=v_mco;
run;

data err4;
  set err3; set mco_est; set eff_marginal4;
  keep m_pred m_l m_v m_w _RMSE_ v v_mco w lambda mxb rmse elsalptmoy alph_m lambda_m elrevptmpoy;
  rmse = _RMSE_;
  elsalptmoy = w*(m_w/m_l);
  alph_m = mxb/rmse;
  lambda_m = PDF('Normal',alph_m)/PROBNORM(alph_m);
  elrevptmpoy = v_mco - v*(lambda/rmse)*lambda_m*(lambda_m+alph_m);
run;

proc print data=err4 noobs;
  var elsalptmoy elrevptmpoy;
  title1 "Elasticité par rapport au salaire";
  title2 "et aux autres revenus";
run;

proc qlim data=chapit19;
     model lfp = enf6 enf618 agef a2 f v chomage /discrete;
     model l = enf6 enf618 agef v w / select(lfp=1);
run;

ods pdf close;