/***************** version révisée du 28/02/2003 *****************/ /*************** pour les valeurs extrêmes **************/ /** Auteur : Nicolas Moreau **/ %macro probit(tab=,y=,x=,id=,points=,mill=,impr1=1,impr2=); data deb; deb=time(); /* pour calcul du temps d'exécution */ cst=1; /* comptabilisation du nombre de périodes par individu */ proc summary data=&tab nway; class &id; output out=nobs(keep=_freq_); data nobs; set nobs; obs=1; data n(keep=n); set nobs; by obs; /* nombre total d'individus */ n+1; if last.obs; run; proc iml; start probit; /************************************************************************/ /***** lecture des données *****/ /************************************************************************/ parm={&x rho}; /* nom des paramètres à estimer */ use &tab; read all var{&y} into y; read all var{&x} into x; close &tab; use n; read all var{n} into n; close n; use nobs; read all var{_freq_} into t; close nobs; tmax=max(t); /* nombre max de périodes */ tmin=min(t); /* nombre min de périodes */ free t; nt=nrow(y); /****************************************************************************/ /**** estimation du Probit sans facteur d'hétérogénéité ****/ /****************************************************************************/ crit=1; lo=0; k=ncol(x); maxiter=2000; bhat=inv(x`*x)*x`*y; /* paramètres initiaux */ do iter=1 to maxiter while(crit>1.0e-6); loglhood=0; z=x*bhat; pdf=exp(-z#z/2)/sqrt(8*atan(1)); /* calcul de la densité f(z) */ cdf=probnorm(z); /* calcul de la fonction de répartition F(z) */ mcdf=probnorm(-z); cond1=z<-5; /* valeurs extrêmes */ cond2=z>5; grad=j(nt,k,0); phi=j(nt,1,0); do i=1 to nt; if (cond1[i,]=1) then do; grad[i,]=(y[i,]#(-z[i,])-(1-y[i,])#(pdf[i,]/(1-cdf[i,])))#x[i,]; /* règle de l'Hôpital */ phi[i,]=-z[i,]#pdf[i,]/(1-cdf[i,]); if (cdf[i,]<=10E-300) then cdf[i,]=10E-300; end; if (cond2[i,]=1) then do; grad[i,]=(y[i,]#(pdf[i,]/cdf[i,])-(1-y[i,])#z[i,])#x[i,]; phi[i,]=z[i,]#pdf[i,]/cdf[i,]; if (mcdf[i,]<=10E-300) then mcdf[i,]=10E-300; end; if (cond1[i,]=0 & cond2[i,]=0) then do; /* calcul du gradient : dérivées premières de la log-vraisemblance/beta */ grad[i,]=(y[i,]#(pdf[i,]/cdf[i,])-(1-y[i,])#(pdf[i,]/(1-cdf[i,])))#x[i,]; phi[i,]=pdf[i,]##2/(cdf[i,]#(1-cdf[i,])); /* voir Davidson et Mackinnon (page 518) */ end; end; grad=grad[+,]`; /* fait la sommation sur toutes les lignes pour chaque colonne et transpose */ inf=(x#phi)`*x; /* matrice d'information, Davidson et Mackinnon (page 519) */ free phi pdf cond1 cond2; loglhood=sum(y#log(cdf)+(1-y)#log(mcdf)); /* normalement, on met log(1-cdf) mais comme distribution symétrique, 1-cdf=mcdf*/ /* comme cela on évite quand cdf proche de 1 des problèmes de calcul */ free cdf mcdf; invinf=inv(inf); free inf; delta=invinf*grad; /* calcul du pas */ /* critères de convergence */ crit1=abs(loglhood-lo); crit2=ssq(grad); crit3=ssq(delta); critall=crit1//crit2//crit3; crit=max(critall); free grad crit1 crit2 crit3 critall; if &impr1=1 then print iter loglhood; /* affichage du nombre d'itérations et de la vraisemblance */ bhat=bhat+delta; /* mise à jour des paramètres estimés */ lo=loglhood; /* mise à jour de la vraisemblance */ end; /* inférence statistique */ etype=sqrt(vecdiag(invinf)); tstat=bhat/etype; pcrit=2#(1-probt(abs(tstat),n-k)); /* impression des résultats */ param=t(parm[,1:k]); print "ESTIMATION DU MODELE PROBIT EMPILE SUR DONNEES DE PANEL",, "NOMBRE D'OBSERVATIONS :" NT, "NOMBRE D'INDIVIDUS :" N, "NOMBRE DE PERIODES MIN et MAX :" tmin tmax, "NOMBRE D'ITERATIONS :" iter, "VALEUR DE LA LOG-VRAISEMBLANCE :" loglhood,, "PARAMETRES ESTIMES", param bhat etype tstat pcrit; /* y réalisés et prédits */ yhat=x*bhat>j(nt,1,0); start predic; y1=sum(y); y0=nt-y1; yhat1=sum(x*bhat[1:k,]>j(nt,1,0)); yhat0=nt-yhat1; y1hat0=sum(y>yhat); y0hat1=sum(y1.0e-6); loglhood=0; inf=0; grad=0; ti_1=0; do i=1 to n; Li=0; gradi=0; ti=t[i,]; ri=2#y[ti_1+1:ti_1+ti,]-j(ti,1,1); /* intégration numérique */ do j=1 to &points; z=ri#(x[ti_1+1:ti_1+ti,]*bhat+teta#abscisse[j,]#j(ti,1,1)); cdf=probnorm(z); pcdf=cdf[#,]; Li=Li+poids[j,]#pcdf; pdf=exp(-z#z/2)/sqrt(8*atan(1)); ratio=j(ti,1,0); /* valeurs extrêmes */ do l=1 to ti; if z[l,]<-5 then ratio[l,]=-z[l,];else ratio[l,]=pdf[l,]/cdf[l,]; end; sgrad=(ri#ratio)`*(x[ti_1+1:ti_1+ti,]||(abscisse[j,]#j(ti,1,1))); free cdf pdf cond ratio; gradi=gradi+poids[j,]#pcdf#sgrad; free pcdf sgrad; end; ti_1=ti_1+ti; grad=grad+(gradi/Li)`; loglhood=loglhood+log(Li/sqrt(4*atan(1))); inf=inf+(gradi/Li)`*(gradi/Li); free gradi Li; end; cov=inv(inf); free inf; delta=cov*grad; /* calcul du pas */ /* critères de convergence */ crit1=abs(loglhood-lo); crit2=ssq(grad); crit3=ssq(delta); critall=crit1//crit2//crit3; crit=max(critall); free grad crit1 crit2 crit3; if &impr2=1 then print iter loglhood; /* affichage du nombre d'itérations et de la vraisemblance */ alpha=alpha+delta; /* mise à jour des paramètres estimés */ bhat=alpha[1:k,]; teta=alpha[k+1,]; lo=loglhood; /* mise à jour de la vraisemblance */ end; /* inférence statistique et impression des résultats */ rho=teta##2/(2+teta##2); /* calcul de l'écart-type de rho par la méthode delta */ etyprho=sqrt(cov[k+1,k+1])#4#teta/(2+teta##2)##2; etype=sqrt(vecdiag(cov[1:k,1:k])); free alpha; etype=etype//etyprho; bhat=bhat//rho; tstat=bhat/etype; pcrit=2#(1-probt(abs(tstat),n-k-1)); param=parm`; print "ESTIMATION DU MODELE PROBIT A EFFETS ALEATOIRES",, "NOMBRE D'OBSERVATIONS :" NT, "NOMBRE D'INDIVIDUS :" N, "NOMBRE DE PERIODES MIN et MAX :" tmin tmax, "NOMBRE D'ITERATIONS :" iter, "VALEUR DE LA LOG-VRAISEMBLANCE :" loglhood,, "PARAMETRES ESTIMES", param bhat etype tstat pcrit; /* y réalisés et prédits */ yhat=x*bhat[1:k,]>j(nt,1,0); run predic; /* calcul et impression des effets marginaux à la moyenne de l'échantillon */ sig2=rho/(1-rho); /* variance de l'effet aléatoire */ moyenne=(j(1,nt,1)*x)/(nt); z=(moyenne*bhat[1:k,])/sqrt(1+sig2); pdfz=exp(-z#z/2)/sqrt(8*atan(1)); effet=pdfz*bhat[1:k,]/sqrt(1+sig2); C=pdfz*i(k)/sqrt(1+sig2)-(bhat[1:k,]/(1+sig2))*z*pdfz*moyenne; var=C*cov[1:k,1:k]*C`; etype=sqrt(vecdiag(var)); tstat=effet/etype; pcrit=2#(1-probt(abs(tstat),n-k-1)); moyenne=moyenne`; param=param[1:k,]; print "ESTIMATION DES EFFETS MARGINAUX A LA MOYENNE DE L'ECHANTILLON", param effet etype tstat pcrit moyenne; free C z pdfz effet etype tstat pcrit moyenne var; /* création de la table incluant les paramètres estimés et de la table*/ /* comprenant la matrice de covariance estimée des paramètres */ bhath=bhat`; create betahath from bhath [colname=parm]; append from bhath; create covarh from cov [colname=parm]; append from cov; free cov etype tstat pcrit crit lo loglhood delta etyprho bhath param; if &mill=1 then do; /*********************************************************************/ /**** calcul des inverses des ratios de Mill *****/ /*********************************************************************/ A1=j(nt,1,0); /* initialisation des vecteurs incluant les ratios */ A2=j(nt,1,0); sig=sqrt(sig2); /* écart-type de l'effet aléatoire */ ti_1=0; do i=1 to n; ti=t[i,]; ri=2#y[ti_1+1:ti_1+ti,]-j(ti,1,1); adenomf=0; espeteps=0; do j=1 to &points; cdf1=probnorm((x[ti_1+1:ti_1+ti,]*bhat[1:k,]+abscisse[j,]#j(ti,1,1))/sig); cdf2=probnorm(ri#(x[ti_1+1:ti_1+ti,]*bhat[1:k,]+abscisse[j,]#j(ti,1,1))/sig); pcdf2=cdf2[#,]; espereta=ri#cdf1/cdf2; free cdf1 cdf2; espeteps=espeteps+poids[j,]#((abscisse[j,]#j(ti,1,1)+espereta)#pcdf2); adenomf=adenomf+poids[j,]#pcdf2; free espereta pcdf2; end; espeteps=espeteps/adenomf; A1[ti_1+1:ti_1+ti,]=(j(1,ti,1)*espeteps/(1+ti#sig##2))#j(ti,1,1); A2[ti_1+1:ti_1+ti,]=espeteps-sig##2#A1[ti_1+1:ti_1+ti,]; ti_1=ti_1+ti; end; free x y; /* création des tables incluant les inverses des ratios de Mill */ /* A1 correspond à celui de l'effet aléatoire et A2 à celui de l'aléa. */ /* Sont ajoutés dans ces tables les identifiants individuels et */ /* temporels */ create A1 from A1; append from A1; create A2 from A2; append from A2; end; finish probit; run probit; quit; data fin; fin=time(); /* pour calcul du temps d'exécution */ cst=1; title "temps d'exécution"; data temps; merge deb fin; by cst; temps=fin-deb; proc print noobs data=temps; var temps; run; %mend probit; /*********************************** exemple d'utilisation **********************************/ /* les données sont dans la table temp, la régressande s'appelle y, les explicatives CST, X1 X2 X3 et X4, l'identificateur */ /* individuel, qui doit être une variable numérique, identn. L'estimation avec 20 points d'évaluation est demandée. On ne */ /* demande ni le calcul des inverses des ratios de Mill généralisés ni l'impression des résultats intermédiaires. */ %probit(tab=temp, y=y, x=cst X1 X2 X3 X4, id=identn, points=20, mill=0, impr1=0, impr2=0);