limited dependent and qualitative variables in econometrics

```{r}

#p.36

library(dummies)

data(iris)

X=as.matrix(iris[,!(colnames(iris) %in% "Species")])

Y=dummy(iris$Species)

m=ncol(Y)

b=array(0.0001,dim=c(ncol(X),m-1))

n=nrow(X)

ite=10^6

eta=0.001

for(l in 1:ite){

mat=cbind(exp(X%*%b),rep(1,n))

P=mat/matrix(rep(apply(mat,1,sum),m),ncol=m)

db=t(Y[,1:(ncol(Y)-1)]-P[,1:(ncol(Y)-1)])%*%X

b=b+eta*t(db)

loglik=sum(Y*log(P))

print(loglik)

}

```


```{r}

#p.46 ordered-response model

#gradient method

library(dummies)

data(iris)

Z=dummy(iris$Species)

X=as.matrix(cbind(rep(1,nrow(iris)),iris[,1:(ncol(iris)-1)]))

class=3

a=c(1:class)/100

b=rep(1,ncol(X))

ite=10^7

eta=10^(-8)

for(l in 1:ite){

b_pre=b;a_pre=a

val=-X%*%b

val=matrix(rep(val,class),ncol=class)

Y=t(t(val)+a)

pthi=dnorm(Y)

bpthi=pnorm(Y)

pthi2=cbind(dnorm(-X%*%b),pthi[,1:(class-1)])

pthi1=pthi

bpthi2=cbind(pnorm(-X%*%b),bpthi[,1:(class-1)])

bpthi1=bpthi

mat=Z*(pthi2-pthi1)/(bpthi1-bpthi2)

db=apply(t(mat)%*%X,2,sum)

da=c()

for(k in 1:class){

da=c(da,sum((Z*pthi1/(bpthi1-bpthi2))[,k]-((Z*pthi2)/(bpthi1-bpthi2))[,k]))

}

b=b+eta*db

a=a+eta*da

a=sort(a)

loglik=sum(Z*log(bpthi1-bpthi2))

#print(sum(abs(b-b_pre))+sum(abs(a-a_pre)))

print(loglik)

}




```

```{r}

#2.16 p.54

library(dplyr)

data=data.frame(num=1:20,y=c(0.20, 0.10, 0.49, 0.26, 0.92, 0.95, 0.18, 0.19, 0.44, 0.79, 0.61, 0.41, 0.49, 0.34, 0.62, 0.99, 0.38, 0.22,0.71,0.40),x1=c(84,87,86,85,82,87,92,94,88,84.9,78,90,88,87,82,84,86,83,85.2,82),x2=c(61,55.5,57,57,50,50,66.5,65,60.5,49.5,49.5,61,59.5,58.4,53.5,54,60,58.8,54,56),x3=c(24.5,24,23,22.5,22,24,25,26,25,24,23,21,26,25.5,24,23,24,24,24,22),sign=0)

for(j in 1:nrow(data)){

data$sign[j]=ifelse(data$y[j]>0.49,1,0)

}

X=as.matrix(cbind(rep(1,nrow(data)),data[,colnames(data) %in% c("x1","x2")]))

Y=data$sign

eta=10^(-0)

beta=rep(0.001,ncol(X))

ite=10^2

P1=0.3;P2=0.4;P3=0.3

for(l in 1:ite){
   
C1=P1/(1-P1-P2)

C2=P2/(1-P1-P2)

R=C1+(1+C1)*exp(X%*%beta)

S=1+C2+C2*exp(X%*%beta)

vec=Y*(1+C1+C2)*exp(X%*%beta)/(R*S)-exp(X%*%beta)/(S*(1+exp(X%*%beta)))
  
dbeta=apply(X*c(vec),2,sum)
 
values=((1+C2)/S+C1/R-1)*Y*(1+C1+C2)*exp(X%*%beta)/(R*S)-((1+C1)/S-exp(X%*%beta)/(1+exp(X%*%beta)))*exp(X%*%beta)/(S*(1+exp(X%*%beta)))

ddbeta=array(0,dim=c(length(beta),length(beta)))

for(i in 1:length(Y)){
   
vector=c(X[i,])   
   
ddbeta=ddbeta+values[i]*t(t(vector))%*%t(vector)   
   
}

beta=beta-eta*solve(ddbeta)%*%dbeta

loglik=sum(Y*log(R)+(1-Y)*log(S)-log(1+exp(X%*%beta)))

print(loglik)
  
}



```

```{r}

#p.98 Morimune(1979)

Breathlessness=matrix(c(9,23,54,121,169,269,404,406,372,7,9,19,48,54,88,117,152,106),ncol=2)
Nobreathlessness=matrix(c(95,105,177,257,273,324,245,225,132,1841,1654,1863,2357,1778,1712,1324,967,526),ncol=2)

#Nobreathlessness[,1]=Nobreathlessness[,1]/10
#Nobreathlessness[,2]=Nobreathlessness[,2]/10



Breath=cbind(Breathlessness,Nobreathlessness)
#Breath=log(Breath)/10

nt=apply(Breath,1,sum)
age_group_und=c(0:8)*5+20
age_group_up=c(0:8)*5+25-1

x=(age_group_und+age_group_up)/2

loglik=function(vector){

a=vector[1];b=vector[2];c=vector[3]

D=1+exp(a*x)+exp(b*x)+exp((a+b+c)*x)

p11=exp((a+b+c)*x)/D

p12=exp((a*x))/D

p21=exp((b*x))/D

p22=1/D

p_mat=cbind(p11,p12,p21,p22)

return(sum(Breath*p_mat))

}


f=function(vector){

a=vector[1];b=vector[2];c=vector[3]

D=1+exp(a*x)+exp(b*x)+exp((a+b+c)*x)

d11a=x*exp((a+b+c)*x)/D-exp((a+b+c)*x)*(x*exp(a*x)+x*exp((a+b+c)*x))/(D^2)
d12a=x*exp(a*x)/D-exp(a*x)*(x*exp(a*x)+x*exp((a+b+c)*x))/(D^2)
d21a=-exp(b*x)*(x*exp(a*x)+x*exp((a+b+c)*x))/(D^2)
d22a=-(x*exp(a*x)+x*exp((a+b+c)*x))/(D^2)

d11b=x*exp((a+b+c)*x)/D-exp((a+b+c)*x)*(x*exp(b*x)+x*exp((a+b+c)*x))/(D^2)
d12b=-exp(a*x)*(x*exp(b*x)+x*exp((a+b+c)*x))/(D^2)
d21b=x*exp(b*x)/D-exp(b*x)*(x*exp(b*x)+x*exp((a+b+c)*x))/(D^2)
d22b=-(x*exp(b*x)+x*exp((a+b+c)*x))/(D^2)

d11c=x*exp((a+b+c)*x)/D-exp((a+b+c)*x)*(x*exp((a+b+c)*x))/(D^2)
d12c=-exp(a*x)*(x*exp((a+b+c)*x))/(D^2)
d21c=-exp(b*x)*(x*exp((a+b+c)*x))/(D^2)
d22c=-(x*exp((a+b+c)*x))/(D^2)

da=sum(Breath*cbind(d11a,d12a,d21a,d22a))
db=sum(Breath*cbind(d11b,d12b,d21b,d22b))
dc=sum(Breath*cbind(d11c,d12c,d21c,d22c))

return(c(da,db,dc))

}

ite=10000

param=rep(10^(-3),3)

eta=10^(-6)

for(l in 1:ite){

vec=param

d=f(param)

param=param+eta*d

print(loglik(param))

}

option=0

if(option==1){

a=param[1];b=param[2];c=param[3]

D=1+exp(a*x)+exp(b*x)+exp((a+b+c)*x)

p11=exp((a+b+c)*x)/D

p12=exp((a*x))/D

p21=exp((b*x))/D

p22=1/D

p=cbind(p11,p12,p21,p22)


}

```





```{r}



#censored and truncated regression models p.149

#The tobit model(censored regression)

#Fair(1977)

loglik=function(beta,sigma){

val1=sum(log(1-pnorm(X[Y==0,]%*%beta,sd=sqrt(sigma))))

val2=sum(Y==1)*log(1/sqrt(2*sigma))

val3=-sum((Y[Y==1]-X[Y==1,]%*%beta)^2)/(2*sigma)

return(val1+val2+val3)

}


data(infert)

colnames(infert)

X=as.matrix(infert[,!(colnames(infert) %in% c("education","case"))])

Y=infert$case

b=rep(0.0001,ncol(X))

ite=10^5

eta=0.0001

for(l in 1:ite){

b_pre=b

sigma=sum((Y[Y==1]-X[Y==1,]%*%b)*Y[Y==1])/length(Y[Y==1])

r0=sqrt(sigma)*dnorm(X[Y==0,]%*%b/sqrt(sigma))/(1-pnorm(X[Y==0,]%*%b/sqrt(sigma)))

b_childer=solve(t(X[Y==1,])%*%X[Y==1,])%*%t(X[Y==1,])%*%Y[Y==1]-sqrt(sigma)*solve(t(X[Y==1,])%*%X[Y==1,])%*%t(X[Y==0,])%*%r0

b=b+eta*(b_childer-b)

print(loglik(b,sigma))

#print(sum((b-b_pre)^2))

}




```

```{r}

#A reparametrization of the tobit model p.156

data(infert)

colnames(infert)

X=as.matrix(infert[,!(colnames(infert) %in% c("education","case"))])

Y=infert$case

ite=30

eta=0.01

b=rep(0.01,ncol(X))

h=0.01

for(l in 1:ite){

db=-t(X[Y==0,])%*%c(dnorm(-X[Y==0,]%*%b)/pnorm(-X[Y==0,]%*%b))+t(X[Y==1,])%*%(h*Y[Y==1]-X[Y==1,]%*%b)

dh=sum(Y==1)/h-sum((h*Y[Y==1]-X[Y==1,]%*%b)*Y[Y==1])

dbb=array(0,dim=c(length(b),length(b)))

X0=X[Y==0,]

for(j in 1:sum(Y==0)){

dbb=dbb+(dnorm(-X[Y==0,]%*%b)/pnorm(-X[Y==0,]%*%b))[1]*(X[Y==0,]%*%b-dnorm(-X[Y==0,]%*%b)/pnorm(-X[Y==0,]%*%b))[1]*t(t(c(X0[j,])))%*%t(c(X0[j,]))

}

X1=X[Y==1,]

dbh=apply(X[Y==1,],2,sum)

for(j in 1:sum(Y==1)){

dbb=dbb-t(t(c(X1[j,])))%*%t(c(X1[j,]))

}

dhh=-sum(Y==1)/(h^2)-sum(Y==1)

dd=rbind(c(dhh,dbh),cbind(dbh,dbb))

d=c(dh,db)

param=c(h,b)

param=param-solve(dd)%*%d

h=param[1];b=param[-1]

print(param)

}

Eyx=X%*%b-sigma*dnorm(d)/pnorm(d)

Vyx=sigma^2-(d*dnorm(d)/pnorm(d)+dnorm(d)^2/(pnorm(d)^2))*sigma^2


```
```{r}

#p.166 Trucated regression model

data("anscombe")

X=as.matrix(anscombe[,colnames(anscombe) %in% c("y1","y2","y3","y4")])

Y=anscombe$x4

L=anscombe$x3

no=ifelse(Y<=L,1,0)

N=length(Y)

loglik=function(b,sig){
   
sig=sqrt(sig)
   
return(-N*log(sqrt(2*pi)*sig)-sum(no*((Y-X%*%b)^2/(sig^2)))/2-sum(log(pnorm((L-X%*%b)/sig))*no))
   
}


#newton-raphson


beta=rep(0.001,ncol(X))

sig2=1

XX=c(sig2,beta)

ite=10^5

eta=10^(-2)


for(l in 1:ite){
   
sig2=XX[1];beta=XX[-1]

print(loglik(beta,sig2))

d=(L-X%*%beta)/sqrt(sig2)
   
dbeta=apply(X*c((Y-X%*%beta)*no),2,sum)/sig2+apply(X*c(no*dnorm(d)/(sqrt(sig2)*pnorm(d))),2,sum)  
   
dsig2=-N/(2*sig2)+sum(no*(Y-X%*%beta)^2)/(2*sig2^2)+sum(no*d*dnorm(d)/pnorm(d))/(2*sig2)

ddsig2=N/(2*sig2^2)-sum(no*(Y-X%*%beta)^2)/(sig2^3)-sum(no*d*dnorm(d)/pnorm(d))/(2*sig2^2)

ddbeta=array(0,dim=c(length(beta),length(beta)))

for(j in 1:nrow(X)){
   
ddbeta=ddbeta-no[j]*t(t(c(X[j,])))%*%t(c(X[j,]))/sig2   
   
}

dsig_beta=-apply(X*c(no*(Y-X%*%beta)),2,sum)/sig2


H=rbind(c(ddsig2,dsig_beta),cbind(dsig_beta,ddbeta))

XX=XX-eta*solve(H)%*%c(dsig2,dbeta)

}

  

```

```{r}

#p.187 Nonnormal error distributions

data("anscombe")

X=as.matrix(anscombe[,colnames(anscombe) %in% c("y1","y2","y3")])

Y=anscombe$y4

alpha=6

loglik=function(beta,sig){
  
mu=-log(alpha)+log(ifelse(X%*%beta>0,X%*%beta,1))-sig^2/2

V=log(Y)-log(ifelse(X%*%beta>0,X%*%beta,1))+sig^2/2

f=pnorm(mu,sd=sig)

val=sum(log(f[Y<alpha]))-sum(log((sig*Y[Y>=alpha])^2))/2-sum(V[Y>=alpha]^2)/2
  
print(val)

}


ite=10^4

eta=0.001

b=rep(0.01,ncol(X))

sigma=10

N1=length(Y[Y<alpha])

N2=length(Y[Y>=alpha])

for(l in 1:ite){
  
mu=-log(alpha)+log(ifelse(X%*%b>0,X%*%b,1))-sigma^2/2  

V=log(Y)-log(ifelse(X%*%b>0,X%*%b,1))+sigma^2/2
  
db1=t(X[Y<alpha,])%*%c((dnorm(mu,sd=sigma)/(pnorm(mu,sd=sigma)*X%*%b))[Y<alpha])  
db2=t(X)%*%c(V/(X%*%b*sigma^2))
  
db=db1+db2

dsigma=-sum((mu*dnorm(mu,sd=sigma)/pnorm(mu,sd=sigma))[Y<alpha])/(sigma^2)-N2/(2*sigma^2)+sum(V[Y>=alpha]^2)/(2*sigma^4)-sum((dnorm(mu,sd=sigma)/pnorm(mu,sd=sigma))[Y<alpha])/2-sum(V[Y>=alpha])/(2*sigma^2)

b=b+eta*db

sigma=sigma+eta*dsigma

print(loglik(b,sigma))

}



```

```{r}

#censored and truncated regression models p.189

#The tobit model(censored regression)

C=1

data(infert)

colnames(infert)

X=as.matrix(infert[,!(colnames(infert) %in% c("education","case"))])

Y=infert$case

loglik=function(beta){

log_val=(X[Y<C,]%*%beta)

log_val=log(ifelse(log_val<0,1,log_val))
  
val=-sum(log(X[Y<C,]%*%beta))-sum(Y[Y<C]/(X[Y<C,]%*%beta)^2)-sum(C/(X[Y==C,]%*%beta)^2)
  
return(val)

}

b=rep(0.001,ncol(X))

ite=10^5

eta=10^(-10)

for(l in 1:ite){

db=-t(X[Y<C,])%*%(1/(X[Y<C,]%*%b))+t(X[Y<C,])%*%(Y[Y<C]/(X[Y<C,]%*%b)^2)+t(X[Y==C,])%*%(C/(X[Y==C,]%*%b)^2)

b=b+eta*db

print(loglik(b))

}


```

```{r}

#p.216

L=c(1000,2000,5000,9000,0,0,4000,0,0,7000)/100

R=c(1,2,5,6,4,8,9,4,3,2)

X1=c(20,45,36,21,84,97,45,23,12,24)

X2=c(69,35,21,14,57,67,1,23,25,45)

n1=length(sum(L>0))

n2=length(sum(L==0))




#Broyden Quasi-Newton Algorithm p.213

f=function(x){

a1=x[1];a2=x[2];b1=x[3];b2=x[4];sig1=x[5];sig2=x[6];p=x[7]  
  
Z1=L-a1*R-b1*X1

Z2=L-a2*R-b2*X2


pnom=ifelse(pnorm(delta)==1,0.9999,pnorm(delta))

delta=((a2-a1)*mean(R)-b1*X1+b2*X2)/(sig1^2+sig2^2-2*p*sig1*sig2)

dsig=sum(((sig1-p*sig2)/(sig1^2+sig2^2-2*p*sig1*sig2))*delta[L==0]*dnorm(delta[L==0])/(1-pnom[L==0]))

dsig1=-n1/sig1-sum(((-2*Z1^3/(sig1^3)+2*p*Z1*Z2/(sig2*sig1^2))/(2*(1-p^2)))[L>0])+dsig

dsig2=-n1/sig2-sum(((-2*Z2^3/(sig2^3)+2*p*Z1*Z2/(sig1*sig2^2))/(2*(1-p^2)))[L>0])+dsig

dp=n1*p/(1-p^2)+(p/((1-p^2)^2))*sum((Z1^2/(sig1^2)+Z2^2/(sig2^2)-2*Z1*Z2/(sig1*sig2))[L>0])+sum((Z1*Z2)[L>0])/((1-p^2)*sig1*sig2)+sum((dnorm(delta)/(1-pnom))[L>0])*(sig1*sig2/(sig1^2+sig2^2-2*p*sig1*sig2))

da1=-n1/(a2-a1)+sum((Z1*R)[L>0]/(sig1^2)-(p*R*Z2)[L>0]/(sig1*sig2))/(1-p^2)-sum((dnorm(delta)/(1-pnom))[L==0])*mean(R)/sqrt(sig1^2+sig2^2-2*p*sig1*sig2)

da2=-n1/(a2-a1)+sum((Z2*R)[L>0]/(sig2^2)-(p*R*Z1)[L>0]/(sig1*sig2))/(1-p^2)-sum((dnorm(delta)/(1-pnom))[L==0])*mean(R)/sqrt(sig1^2+sig2^2-2*p*sig1*sig2)

db1=sum(((Z1/(sig1^2)-p*Z2/(sig1*sig2))*X1)[L>0])/(1-p^2)-sum(((dnorm(delta)/(1-pnom))*X1)[L==0])/(sig1^2+sig2^2-2*p*sig1*sig2)

db2=sum(((Z2/(sig2^2)-p*Z1/(sig1*sig2))*X2)[L>0])/(1-p^2)-sum(((dnorm(delta)/(1-pnom))*X2)[L==0])/(sig1^2+sig2^2-2*p*sig1*sig2)

return(c(da1,da2,db1,db2,dsig1,dsig2,dp))

}

X=c(sample(c(1:9),6),0.5)

ite=10000

eta=10^(-4)

for(l in 1:ite){
  
X_pre=X  
  
X=X-eta*f(X)  

#print(sum(abs(X-X_pre)))

a1=X[1];a2=X[2];b1=X[3];b2=X[4];sig1=X[5];sig2=X[6];p=X[7]  

delta=((a2-a1)*mean(R)-b1*X1+b2*X2)/(sig1^2+sig2^2-2*p*sig1*sig2)

loglik=-n1*log(sig1)-n1*log(sig2)-n1*log(1-p^2)/2-sum((Z1^2/(sig1^2)+Z2^2/(sig2^2)-2*Z1*Z2/(sig1*sig2))[L>0])/(2*(1-p^2))+sum(log(1-ifelse(pnorm(delta)==1,0.9999,pnorm(delta)))[L==0])

print(loglik)
  
}



```