R: maximum entropy regression linear and non linear

29Dec11

This code is for implementation of maximum entropy  non linear and linear estimation of regression in R. You will need Alabama package.

## PARA REGRESIONES LINEALES.
## Funcion GENERALIZED MAXIMUM ENTROPY

gme.v1=function(Y=Y,X=X,M=M,C=C,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
  ## Support Matrix of dim k x M.

if (opt2==”noprior”) {   
    z1=rep(0,M)
    nz1=floor(length(z1)/2)
    dum<-1/nz1
    dum<-seq(0,1,by=(dum))*C
    z1<-c(-rev(dum[2:(nz1+1)]),dum)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## inicial po.
    p1o=rep(1,length(z1))
    p1o=p1o/sum(p1o)
    po<-rep(p1o,dim(X)[2])

}

if (opt2==”prior”) {
    z1=rep(0,M)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## rownames
    dum1<-rep(“beta”,dim(X)[2])
    dum2<-seq(1:dim(X)[2])
    dum1<-paste(dum1,dum2,sep=”")
    rownames(Z)<-dum1
    ## colnames
    dum1<-rep(“p”,M)
    dum2<-seq(1:M)
    dum1<-paste(dum1,dum2,sep=”")
    dum2<-rep(dum1,dim(X)[2])
    colnames(Z)<-dum2
    ##
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
        dum=dum+(M)
    }
   
   
    ## inicial po.
    p1o=rep(0,length(z1))
    po<-rep(p1o,dim(X)[2])
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        po[dum1:(dum1+(M-1))]<-matprior2[i,]
        dum=dum+M
    }

}
   
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)

V=kronecker(diag(1,length(Y)), t(v1))

## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))

## Valor inicial total.
pwo<-c(po,wo)

## k
k=length(po)

## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    ## tama?o de restricciones.
    res.1<-dim(V)[2]
    ## restricciones de probabilidades.
    res.2<-2
    res<-res.1+res.2
    h<-rep(NA,res)
    ## Restricciones lineales.
    h[1:res.1]<-X%*%(Z%*%p)+(V%*%w)-Y
    ## Restricciones de probabilidades.
    h[(res.1+1)]<-sum(p,na.rm=TRUE)-1
    h[(res.1+2)]<-sum(w,na.rm=TRUE)-1
    h
}   

## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    h<-rep(NA,(2*len))
    ## Restricciones lineales.
    h[1:len]<-pw
    h[(len+1):(2*len)]<- -pw+1
}   

H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
    -Objective
}

## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Grad.p<-log(p)+1
    Grad.w<-log(w)+1
    Grad<-c(Grad.p,Grad.w)
    Grad
}
## El default de auglag y optim es minimizar.

################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
    results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}

if (opt1==”gr”) {
    ## con gradiente
    results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]

e<-V%*%results$par[(k+1):length(results$par)]

    list(beta=beta,e=e)
}

## Generalized Maximum Entropy no lineal  Especializada para

gme.nl=function(Y=Y,X=X,M=M,C=C,fun=fun,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
## inicializacion de I y L
  I=X[,"I"]
  L=X[,"L"]
 
## Support Matrix of dim k x M.

if (opt2==”noprior”) { 
    z1=rep(0,M)
    nz1=floor(length(z1)/2)
    dum<-1/nz1
    dum<-seq(0,1,by=(dum))*C
    z1<-c(-rev(dum[2:(nz1+1)]),dum)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## inicial po.
    p1o=rep(1,length(z1))
    p1o=p1o/sum(p1o)
    po<-rep(p1o,dim(X)[2])

}

if (opt2==”prior”) {
    z1=rep(0,M)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## rownames
    dum1<-rep(“beta”,dim(X)[2])
    dum2<-seq(1:dim(X)[2])
    dum1<-paste(dum1,dum2,sep=”")
    rownames(Z)<-dum1
    ## colnames
    dum1<-rep(“p”,M)
    dum2<-seq(1:M)
    dum1<-paste(dum1,dum2,sep=”")
    dum2<-rep(dum1,dim(X)[2])
    colnames(Z)<-dum2
    ##
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
        dum=dum+(M)
    }
   
   
    ## inicial po.
    p1o=rep(0,length(z1))
    po<-rep(p1o,dim(X)[2])
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        po[dum1:(dum1+(M-1))]<-matprior2[i,]
        dum=dum+M
    }

}
   
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)

V=kronecker(diag(1,length(Y)), t(v1))

## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))

## Valor inicial total.
pwo<-c(po,wo)

## k
k=length(po)

## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
  len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    ## tama?o de restricciones.
    res.1<-dim(V)[2]
    ## restricciones de probabilidades.
    res.2<-2
    res<-res.1+res.2
    h<-rep(NA,res)
    ## Restricciones lineales.
    theta<-(Z%*%p)
    h[1:res.1]<-fun(theta=theta,X=X)+(V%*%w)-Y
    ## Restricciones de probabilidades.
    h[(res.1+1)]<-sum(p,na.rm=TRUE)-1
    h[(res.1+2)]<-sum(w,na.rm=TRUE)-1
    h
}   

## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    h<-rep(NA,(2*len))
    ## Restricciones lineales.
    h[1:len]<-pw
    h[(len+1):(2*len)]<- -pw+1
}   

H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
    -Objective
}

## El default de auglag y optim es minimizar.
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Grad.p<-log(p)+1
    Grad.w<-log(w)+1
    Grad<-c(Grad.p,Grad.w)
    Grad
}

## El default de auglag y optim es minimizar.

################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
    results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}

if (opt1==”gr”) {
    ## con gradiente
    results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]

e<-V%*%results$par[(k+1):length(results$par)]

    list(beta=beta,e=e)
}

#gme.nl=cmpfun(gme.nl)
##
##
## Generalized Maximum Entropy no lineal v2 especializado para stock de capital.

gme.nl2=function(Y=Y,X=X,M=M,C=C,fun=fun,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
## inicializacion de I y L
  I=X[,"I"]
  L=X[,"L"]
 
## Support Matrix of dim k x M.

if (opt2==”noprior”) { 
  z1=rep(0,M)
    nz1=floor(length(z1)/2)
    dum<-1/nz1
    dum<-seq(0,1,by=(dum))*C
    z1<-c(-rev(dum[2:(nz1+1)]),dum)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## inicial po.
    p1o=rep(1,length(z1))
    p1o=p1o/sum(p1o)
    po<-rep(p1o,dim(X)[2])

}

if (opt2==”prior”) {
    z1=rep(0,M)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## rownames
    dum1<-rep(“beta”,dim(X)[2])
    dum2<-seq(1:dim(X)[2])
    dum1<-paste(dum1,dum2,sep=”")
    rownames(Z)<-dum1
    ## colnames
    dum1<-rep(“p”,M)
    dum2<-seq(1:M)
    dum1<-paste(dum1,dum2,sep=”")
    dum2<-rep(dum1,dim(X)[2])
    colnames(Z)<-dum2
    ##
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
        dum=dum+(M)
    }
   
   
    ## inicial po.
    p1o=rep(0,length(z1))
    po<-rep(p1o,dim(X)[2])
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        po[dum1:(dum1+(M-1))]<-matprior2[i,]
        dum=dum+M
    }

}

 
## Matriz para obtener restricciones de probabilidades
z.id=rep(1,M)
Idp=kronecker(diag(1,dim(X)[2]), t(z.id))
   
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)

V=kronecker(diag(1,length(Y)), t(v1))

 
## Matriz para obtener restricciones de probabilidades
v1.id=rep(1,length(v1))
Idw=kronecker(diag(1,length(Y)), t(v1.id))
 
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))

## Valor inicial total.
pwo<-c(po,wo)

## k
k=length(po)

## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
  len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    ## tama?o de restricciones.
    res.1<-length(Y)
    ## restricciones de probabilidades.
    res.2<-dim(Z)[1]
    res<-res.1+res.2+res.1
    h<-rep(NA,res)
    ## Restricciones lineales.
    theta<-(Z%*%p)
    h[1:res.1]<-fun(theta=theta,X=X)+(V%*%w)-Y
    ## Restricciones de probabilidades.
    h[(res.1+1):(res.1+(res.2))]<-as.numeric(Idp%*%p-rep(1,times=(res.2/2)))
    h[(res.1+(res.2)+1):length(h)]<-as.numeric(Idw%*%w-rep(1,times=res.1))
    h
}   

## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    h<-rep(NA,(2*len))
    ## Restricciones lineales.
    h[1:len]<-pw
    h[(len+1):(2*len)]<- -pw+1
}   

H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
    -Objective
}

## El default de auglag y optim es minimizar.
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Grad.p<-log(p)+1
    Grad.w<-log(w)+1
    Grad<-c(Grad.p,Grad.w)
    Grad
}

## El default de auglag y optim es minimizar.

################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
    results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}

if (opt1==”gr”) {
    ## con gradiente
    results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]

e<-V%*%results$par[(k+1):length(results$par)]

Hn=-results$value/log(length(results$par))

  list(beta=beta,e=e,Hn=Hn)
}

#gme.nl=cmpfun(gme.nl)

gme.v2=function(Y=Y,X=X,M=M,C=C,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2) {
  ## Support Matrix of dim k x M.

if (opt2==”noprior”) { 
    z1=rep(0,M)
    nz1=floor(length(z1)/2)
    dum<-1/nz1
    dum<-seq(0,1,by=(dum))*C
    z1<-c(-rev(dum[2:(nz1+1)]),dum)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## inicial po.
    p1o=rep(1,length(z1))
    p1o=p1o/sum(p1o)
    po<-rep(p1o,dim(X)[2])

}

if (opt2==”prior”) {
    z1=rep(0,M)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## rownames
    dum1<-rep(“beta”,dim(X)[2])
    dum2<-seq(1:dim(X)[2])
    dum1<-paste(dum1,dum2,sep=”")
    rownames(Z)<-dum1
    ## colnames
    dum1<-rep(“p”,M)
    dum2<-seq(1:M)
    dum1<-paste(dum1,dum2,sep=”")
    dum2<-rep(dum1,dim(X)[2])
    colnames(Z)<-dum2
    ##
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
        dum=dum+(M)
    }
   
   
    ## inicial po.
    p1o=rep(0,length(z1))
    po<-rep(p1o,dim(X)[2])
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        po[dum1:(dum1+(M-1))]<-matprior2[i,]
        dum=dum+M
    }

}
   
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)

V=kronecker(diag(1,length(Y)), t(v1))

## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))

## Valor inicial total.
pwo<-c(po,wo)

## k
k=length(po)
 
  ##  matrices de restricciones
  ## Matriz para obtener restricciones de probabilidades
  v1.id=rep(1,length(v1))
  Idw=kronecker(diag(1,length(Y)), t(v1.id))
  ## Matriz para obtener restricciones de probabilidades
  z.id=rep(1,M)
  Idp=kronecker(diag(1,dim(X)[2]), t(z.id))

## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    ## tama?o de restricciones.
    res.1<-length(Y)
    ## restricciones de probabilidades.
    res.2<-dim(Z)[1]
    res<-res.1+res.2+res.1
    h<-rep(NA,res)
    ## Restricciones lineales.
    h[1:res.1]<-X%*%(Z%*%p)+(V%*%w)-Y
    ## Restricciones de probabilidades.
    ## Restricciones de probabilidades.
  h[(res.1+1):(res.1+(res.2))]<-as.numeric(Idp%*%p-rep(1,times=(res.2/2)))
    h[(res.1+(res.2)+1):length(h)]<-as.numeric(Idw%*%w-rep(1,times=res.1))
    h
}   

## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    h<-rep(NA,(2*len))
    ## Restricciones lineales.
    h[1:len]<-pw
    h[(len+1):(2*len)]<- -pw+1
}   

H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
    -Objective
}

## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Grad.p<-log(p)+1
    Grad.w<-log(w)+1
    Grad<-c(Grad.p,Grad.w)
    Grad
}
## El default de auglag y optim es minimizar.

################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
    results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}

if (opt1==”gr”) {
    ## con gradiente
    results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]

e<-V%*%results$par[(k+1):length(results$par)]

Hn<–results$value

    list(beta=beta,e=e,Hn=Hn)
}

#———————————————————————-
# entropia normalizada, recordad que la distribucion uniforme maximiza la entropia 1/n, entonces
# se normaliza la entropia dividiendo para ln n . Hn = H/ln n
# -results$fn/log(dim(X)[1])
## Generalized Maximum Entropy no lineal v2 especializado para stock de capital.

gme.nl3=function(Y=Y,X=X,M=M,C=C,fun=fun,opt1=opt1,opt2=opt2,matprior=matprior,matprior2=matprior2,fix.t=fix.t) {
## inicializacion de I y L
  I=X[,"I"]
  L=X[,"L"]
 
## Support Matrix of dim k x M.

if (opt2==”noprior”) { 
  z1=rep(0,M)
  nz1=floor(length(z1)/2)
    dum<-1/nz1
    dum<-seq(0,1,by=(dum))*C
    z1<-c(-rev(dum[2:(nz1+1)]),dum)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## inicial po.
    p1o=rep(1,length(z1))
    p1o=p1o/sum(p1o)
    po<-rep(p1o,dim(X)[2])

}

if (opt2==”prior”) {
    z1=rep(0,M)
    Z=kronecker(diag(1,dim(X)[2]), t(z1))
    ## rownames
    dum1<-rep(“beta”,dim(X)[2])
    dum2<-seq(1:dim(X)[2])
    dum1<-paste(dum1,dum2,sep=”")
    rownames(Z)<-dum1
    ## colnames
    dum1<-rep(“p”,M)
    dum2<-seq(1:M)
    dum1<-paste(dum1,dum2,sep=”")
    dum2<-rep(dum1,dim(X)[2])
    colnames(Z)<-dum2
    ##
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        Z[i,dum1:(dum1+(M-1))]<-matprior[i,]
        dum=dum+(M)
    }
   
   
    ## inicial po.
    p1o=rep(0,length(z1))
    po<-rep(p1o,dim(X)[2])
    dum=0
    for (i in 1:dim(X)[2]) {
        dum1=1+dum
        po[dum1:(dum1+(M-1))]<-matprior2[i,]
        dum=dum+M
    }

}

 
## Matriz para obtener restricciones de probabilidades
z.id=rep(1,M)
Idp=kronecker(diag(1,dim(X)[2]), t(z.id))
   
## Soportes para w, en este caso se va a utilizar la
## regla 3sigma.
## sigma es la desviaci?n est?ndar emp?rica de Y.
sigma=sd(Y)
## El soporte para J=5.
## El valor m?nimo debe ser -3*sigma
## El valor m?ximo debe ser 3*sigma
## Deben ser sim?tros alrrededor de cero.
## vt=(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)
v1<-c(-3*sigma,-1.5*sigma,0,1.5*sigma,3*sigma)

V=kronecker(diag(1,length(Y)), t(v1))

 
## Matriz para obtener restricciones de probabilidades
v1.id=rep(1,length(v1))
Idw=kronecker(diag(1,length(Y)), t(v1.id))
 
## inicial po.
w1o=rep(1,length(v1))
w1o=w1o/sum(w1o)
wo<-rep(w1o,length(Y))

## Valor inicial total.
pwo<-c(po,wo)

## k
k=length(po)

## Equality restrictions.
## h[i]=0
heq<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
  len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    ## tama?o de restricciones.
    res.1<-length(Y)
    ## restricciones de probabilidades.
    res.2<-dim(Z)[1]
    res<-res.1+res.2+res.1
    h<-rep(NA,res)
    ## Restricciones lineales.
    theta<-(Z%*%p)
    h[1:res.1]<-fun(theta=theta,X=X)+(V%*%w)-Y
    ## Restricciones de probabilidades.
    h[(res.1+1):(res.1+(res.2))]<-as.numeric(Idp%*%p-rep(1,times=(res.2/2)))
    h[(res.1+(res.2)+1):length(h)]<-as.numeric(Idw%*%w-rep(1,times=res.1))
   
  # fix.t
  hf<-as.numeric(V%*%w)[fix.t]
  h<-c(h,hf)
 
  h
}   

## Inequality restrictions.
## h[i]>0
hin<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    h<-rep(NA,(2*len))
    ## Restricciones lineales.
    h[1:len]<-pw
    h[(len+1):(2*len)]<- -pw+1
}   

H.obj<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Objective=as.numeric(-(t(p)%*%log(p)))+as.numeric(-(t(w)%*%log(w)))
    -Objective
}

## El default de auglag y optim es minimizar.
## gradiente, recibe R^n y bota R^n.
H.gr<-function(pw=pw,Z=Z,V=V,k=k,Y=Y,fun=fun,X=X) {
## pw vector de par?metros que concatena a p y w.
## Z, es la matriz de soportes de los Beta.
## V, es la matriz de soportes de los errores.
## k, longitud del par?metro p.
    len=length(pw)
    p=pw[1:k]
    w=pw[(k+1):len]
    Grad.p<-log(p)+1
    Grad.w<-log(w)+1
    Grad<-c(Grad.p,Grad.w)
    Grad
}

## El default de auglag y optim es minimizar.

################################################################
## Optimizaci?n.
################################################################
if (opt1==”nongr”) {
    results<-constrOptim.nl(par=pwo, fn=H.obj, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}

if (opt1==”gr”) {
    ## con gradiente
    results<-constrOptim.nl(par=pwo, fn=H.obj, gr=H.gr, heq=heq,hin=hin,Z=Z,V=V,k=k,Y=Y,X=X,fun=fun)
}
##beta
## recobro los p del par=pw y multiplico por el soporte Z, para obtener el vector de coeficientes.
beta<-Z%*%results$par[1:k]

e<-V%*%results$par[(k+1):length(results$par)]

Hn=-results$value/log(length(results$par))

  list(beta=beta,e=e,Hn=Hn)
}

Advertisement

Filed under: R   |  Leave a Comment