/*
This program estimates a simultaneous sequential tobit model
in reduced form. The model is set up as follows.
1. Probit-part: This is the selection equation. The Tobit-part
is only observed if the binary variable is one.
2. Tobit-part: This is the levels equation.
The program runs a Full Information Maximum Likelihood procedure.
Recommendation: Start iterations in the full model with BFGS and
sitch then to BHHH.
REFERENCE:
Kaiser, U. (2001): Product Innovation and Product Innovation
Marketing: Theory and Microeconometric Evidence. ZEW mimeo.
*/
/* prerequisites */
new;
library maxlik, pgraph;
#include maxlik.ext;
#include indices.src;
maxset;
format /rd 12,4 ;
/* set global MAXLIK variables */
_max_Algorithm = 5 ; /* see MAXLIK manual, p. 42 */
_max_LineSearch = 1 ; /* see MAXLIK manual, p. 45 */
_max_CovPar = 2 ; /* see MAXLIK manual, p. 42 */
_max_gradtol = 1e-5 ; /* see MAXLIK manual, p. 43 */
_max_GradMethod = 1 ; /* see MAXLIK manual, p. 43 */
_max_MaxIters = 500 ; /* see MAXLIK manual, p. 45 */
output file = c:\fabel\progs\simult.out reset ;
/* COMPLETE DATA SET */
open f1 = c:\fabel\data\mdp.dat varindxi ;
mat = readr(f1, rowsf(f1));
/* all variables */
mat = mat[.,ipd ]~mat[.,iiama_ia]~
mat[.,iH ]~mat[.,idlexe ]~
mat[.,ihandel ]~mat[.,iverk ]~
mat[.,isonst ]~
mat[.,iln_l ]~mat[.,iln_l2 ]~
mat[.,iost ]~
mat[.,im_koop ]~
ones(rows(mat),1)~mat[.,is_high ]~mat[.,idive_br];
names = { pd , iama_ia,
H , dlexe,
handel , verk,
sonst ,
ln_l , lnl2,
ost ,
m_koop ,
cons ,
s_high , dive_br };
exclsel = 1;
excl = 2;
/* FULL DATA SET */
create full = c:\fabel\data\full with ^names,0,8;
call writer(full,mat);
/* ==================================== */
/* STARTING VALUES: BINARY PROBIT MODEL */
/* ==================================== */
namesn = names[1]|names[3:rows(names)];
_max_ParNames = namesn[2:rows(namesn)];
b0 = zeros(rows(_max_ParNames),1);
{binstart,f,g,cov,ret}=maxlik("c:\\fabel\\data\\full",namesn,&probit,b0);
call maxprt(binstart,f,g,cov,ret) ;
binnam = _max_ParNames;
/* hier
/* ==================================== */
/* STARTING VALUES: TYPE II TOBIT */
/* ==================================== */
x = mat[.,3:cols(mat)];
lambda = pdfn(x*binstart)./cdfn(x*binstart);
x = x~lambda~mat[.,2];
x = selif(x,x[.,cols(x)].>0);
y = x[.,cols(x)];
x = x[.,1:cols(x)-1-excl];
blevel = inv(x'*x)*x'y;
blevel = blevel[1:rows(blevel)-1];
b0 = binstart|blevel|0.73|-.11;
s2 = { s2 };
s12 = { s12 };
_max_ParNames = names[3:rows(names)]|names[3:rows(names)-excl]|s2|s12;
_max_CovPar = 2 ; /* see MAXLIK manual, p. 42 */
_max_MaxIters = 10000 ; /* see MAXLIK manual, p. 45 */
_max_Algorithm = 2 ; /* see MAXLIK manual, p. 42 */
_max_GradTol = 1e-5 ; /* see MAXLIK manual, p. 42 */
{bsim,f1,g1,cov1,ret1}=maxlik("c:\\rigorosu\\fabel\\data\\full",names,&tobit2,b0);
call maxprt(bsim,f1,g1,cov1,ret1) ;
hv1 = (rows(bsim)+excl-2)/2;
bsel = bsim[1:hv1];
blev = bsim[hv1+1:hv1+hv1-excl];
covsel = cov1[1:hv1,1:hv1];
covlev = cov1[hv1+1:hv1+hv1-excl,hv1+1:hv1+hv1-excl];
sqrt(diag(covlev));
format /rd 12,4 ;
/* ================================================ */
print"/* Wald-tests for joint significance */";
/* ================================================ */
print"Firm size";
s1 = 6; s2 = 7;
print"Selection equation";
theta = bsel[s1:s2];
var = covsel[s1:s2,s1:s2];
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
print"Level equation";
theta = blev[s1:s2];
var = covlev[s1:s2,s1:s2];
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
print"Sector dummies";
s1 = 3; s2 = 5;
print"Selection equation";
theta = bsel[s1:s2];
var = covsel[s1:s2,s1:s2];
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
print"Level equation";
theta = blev[s1:s2];
var = covlev[s1:s2,s1:s2];
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
print"Entire specification";
print"Selection equation";
theta = bsel[1:rows(bsel)-excl-1]|bsel[rows(bsel)-1:rows(bsel)];
var =
(covsel[1:rows(bsel)-excl-1,1:rows(bsel)-excl-1]~covsel[1:rows(bsel)-excl-1,rows(bsel)-1:rows(bsel)])|
(covsel[rows(bsel)-1:rows(bsel),1:rows(bsel)-excl-1]~covsel[rows(bsel)-1:rows(bsel),rows(bsel)-1:rows(bsel)]);
sqrt(diag(var ));
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
print"Level equation";
theta = blev[1:rows(blev)-1];
var = covlev[1:rows(blev)-1,1:rows(blev)-1];
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
print"Entire specification";
theta = bsel[1:rows(bsel)-excl-1]|bsel[rows(bsel)-excl+1:rows(bsel)]|
blev[1:rows(blev)-1];
var = (cov1[1:9,1:9]~cov1[1:9,11:12]~cov1[1:9,13:21])|
(cov1[11:12,1:9]~cov1[11:12,11:12]~cov1[11:12,13:21])|
(cov1[13:21,1:9]~cov1[13:21,11:12]~cov1[13:21,13:21]);
R = eye(rows(theta));
W = (R*theta)'*inv(R*var*R')*(R*theta);
W~cdfchic(W,rows(R));
/* ================================================ */
print"/* Two-step results */";
/* ================================================ */
_max_Algorithm = 5 ; /* see MAXLIK manual, p. 42 */
namn = names[1]|names[3:rows(names)-excl-exclsel-1]|names[rows(names)-excl:rows(names)];
namn = namn[1:7]|namn[9:rows(namn)];
matn = mat[.,1]~mat[.,3:cols(mat)-excl-exclsel-1]~mat[.,cols(mat)-excl:cols(mat)];
matn = matn[.,1:7]~matn[.,9:cols(matn)];
mhat = mat[.,3:cols(mat)-excl]*blev;
matn = matn~mhat;
hv = { mhat };
namn = namn|hv ;
closeall;
create fulln = c:\rigorosu\fabel\data\fulln with ^namn,0,8;
call writer(fulln,matn);
namesn = namn;
_max_ParNames = namesn[2:rows(namesn)];
b0 = zeros(rows(_max_ParNames),1)+.1;
{binstart,f,g,cov,ret}=maxlik("c:\\rigorosu\\fabel\\data\\fulln",namesn,&probit,b0);
call maxprt(binstart,f,g,cov,ret) ;
hier */
/* SIMULTANEOUS MODEL */
proc tobit2(b,dat);
local p,e,u,hv1,hv2,x,M,D0,D10,D1M,pipd,pim,thetax,thetaM,sigma,rho,dp,blau,gelb,rosa,gruen;
x = dat[.,1];
M = dat[.,2];
D0 = x.==0;
D10 = (x.==1) .and (M.==0);
D1M = (x.==1) .and (M.>0);
pipd = dat[.,3:cols(dat)];
pim = dat[.,3:cols(dat)-excl];
hv1 = (rows(b)+excl-2)/2;
thetax = b[1:hv1];
thetaM = b[hv1+1:hv1+hv1-excl];
sigma = b[rows(b)-1];
rho = b[rows(b)];
blau = ln(cdfn(-pipd*thetax));
gelb = ln(cdfbvn((-pim*thetaM)/sigma,pipd*thetax,-rho));
gruen = ln(cdfn( ( pipd*thetax+(M-pim*thetam)*rho/sigma )/sqrt(1-rho^2) ) )
-.5*((M-pim*thetaM)/sigma)^2-ln(sqrt(2*Pi)*sigma);
retp( ( D0.*blau + D10.*gelb + D1M.*gruen ) );
endp;
/* BINARY PROBIT PROCEDURE AND GRADIENT */
proc probit(b,dat);
local y,x,p,dp;
y = dat[.,1];
x = dat[.,2:cols(dat)];
p = cdfn(x*b);
dp = p.>1e-18;
p = p.*dp+(1-dp)*1e-18;
retp( y.*ln(p) + (1-y).*ln(1-p) );
endp;
proc probitg(b,dat);
local y,x,p,dp;
y = dat[.,1];
x = dat[.,2:cols(dat)];
p = cdfn(x*b);
dp = p.>1e-18;
p = p.*dp+(1-dp)*1e-18;
retp( (y-p)./(p.*(1-p)).*pdfn(x*b).*x );
endp;
output off;