Ricardian Constant Cost Industry

# Problem Set 9.a -----------------------------------------------------------

### Mixed Complementarity Problem [MCP]

#system of equations to solve

rm(list=ls()) #will remove ALL objects



fcnuMax = function(inputs){
  
  PX=1

  XDH=inputs[1]
  YDH=inputs[2]
  XDF=inputs[3]
  YDF=inputs[4]
  XSH=inputs[5]
  YSH=inputs[6]
  XSF=inputs[7]
  YSF=inputs[8]
  PY=inputs[9]
  lam1=inputs[10]
  lam2=inputs[11]
  lam3=inputs[12]
  lam4=inputs[13]

  
  # Utility  -------------------------------------------
  
            
  
  res1=XDH*((ALPHA)*(XDH^(ALPHA-1))*(YDH^(BETA))-PX*lam3)
  res2=YDH*((BETA)*(XDH^(ALPHA))*(YDH^(BETA-1))-PY*lam3)        
  res3=XDF*((GAMA)*(XDF^(GAMA-1))*(YDF^(TETA))-PX*lam4)  
  res4=YDF*((TETA)*(XDF^(GAMA))*(YDF^(TETA-1))-PY*lam4) 
  
  

  # Production  -------------------------------------------
  
  
   res5=XSH*(5*lam1+lam3*PX)
   res6=YSH*(-lam1+lam3*PY)
   res7=XSF*(-(1/5)*lam2+lam4*PX)
   res8=YSF*(-lam2+lam4*PY)
   res9=lam1*(YSH-500+5*XSH)
   res10=lam2*(YSF-100+(1/5)*XSF)
   res11=lam3*(PX*XSH+PY*YSH-PX*XDH-PY*YDH)
   res12=lam4*(PX*XSF+PY*YSF-PX*XDF-PY*YDF)
   
   
  
  # Market Clearing Condition
  
 
  res13=PY*(YSH+YSF-YDH-YDF)
  res14=PX*(XSH+XSF-XDH-XDF) # ignore this


  return(c(res1,res2,res3,res4,res5,res6,res7,res8,res9,res10,res11,res12,res13)) 
  
}


# Parameters
PX=1
ALPHA=1/3
BETA=2/3
GAMA=1/3
TETA=2/3


library("alabama") #load "alabama" Package

## Loading required package: numDeriv

library("rootSolve") #load "rootSolve" Package

## 
## Attaching package: 'rootSolve'

## The following object is masked from 'package:numDeriv':
## 
##     hessian

library("numDeriv") #load "numDeriv" Package

# XDH YDH XDF YDF XSH YSH XSF YSF py lam1 lam2 lam3 lam4
stval=c(333.33, 333.33, 166.67, 166.67, 0, 500, 500, 0, 2, 0.67, 1.67, 0.33, 0.33)

solUFOC=multiroot(fcnuMax, stval)
solinput=solUFOC$root
startF = paste(round(solinput,2), collapse=", ")
startF

## [1] "333.33, 333.33, 166.67, 166.67, 0, 500, 500, 0, 2, 0.67, 1.67, 0.33, 0.33"

PX=1
YDHstar=solinput[1];YDHstar

## [1] 333.3333

XDHstar=solinput[2];XDHstar

## [1] 333.3333

YDFstar=solinput[3];YDFstar

## [1] 166.6667

XDFstar=solinput[4];XDFstar

## [1] 166.6667

YSHstar=solinput[5];YSHstar

## [1] 0

XSHstar=solinput[6];XSHstar

## [1] 500

YSFstar=solinput[7];YSFstar

## [1] 500

XSFstar=solinput[8];XSFstar

## [1] -8.142197e-22

PYstar=solinput[9];PYstar

## [1] 2