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


# ******** Analytically solving for autarky and free traded equilibriums ****

# Parameters 

A0U = 25     # supply intercept country U
A1U = 1.25   # supply slope country U
B0U = 250    # demand intercept country U
B1U = -2.5   # demand slope country U

A0E = 50    # supply intercept country E
A1E = 2     # supply slope country E
B0E = 125   # demand intercept country E
B1E = -1.75 # demand slope country E


PU_Au = (-A0U + B0U)/(A1U - B1U); PU_Au # Autarky price for country U
## [1] 60
PE_Au = (-A0E + B0E)/(A1E - B1E); PE_Au # Autarky price for country E
## [1] 20
QU_Au = (-A0U*B1U + A1U*B0U)/(A1U - B1U); QU_Au # Autarky Quantity for country U
## [1] 100
QE_Au = (-A0E*B1E + A1E*B0E)/(A1E - B1E); QE_Au # Autarky Quantity for country E
## [1] 90
PW = ((A0U - B0U) + (A0E - B0E))/((-A1U + B1U) + (-A1E + B1E)); PW # World Price
## [1] 40
QE_ES = ((A0E - B0E - B0U + A0U)*(A1E - B1E) + (A0E - B0E)*
         (-A1U + B1U - A1E + B1E))/((-A1U + B1U) + (-A1E + B1E)); QE_ES # Excess supply of country U (Note ES=ED)
## [1] 75
QU_ED =  QE_ES # Excess demand of country E (Note ES=ED)

# ******** Numerically solving for autarky and free traded equilibriums ****
 
fcnSDU = function(x){
  QU = x[1]
  PU = x[2]
  
  SupplyU = QU - (A0U + A1U*PU) # QU - (A0U + A1U*PU) = 0 in equilibrium
  DemandU = QU - (B0U + B1U*PU) # QU - (B0U + B1U*PU) = 0 in equilibrium
  return(c(SupplyU,DemandU))
}

fcnSDE = function(x){
  QE = x[1]
  PE = x[2]

  SupplyE = QE - (A0E + A1E*PE) # QE - (A0E + A1E*PE) = 0 in equilibrium
  DemandE = QE - (B0E + B1E*PE) # QE - (B0E + B1E*PE) = 0 in equilibrium
  return(c(SupplyE,DemandE))
}

fcnESEDw = function(x){
  QW = x[1]
  PW = x[2]

  ESupplyw = QW - ((B0U + B1U*PW) - (A0U + A1U*PW))
  EDemandw = QW - ((A0E + A1E*PW) - (B0E + B1E*PW))
  
  return(c(ESupplyw,EDemandw))
}

library("rootSolve")
solnU = multiroot(fcnSDU, c(1, 1))
QU = solnU$root[1]; QU
## [1] 99.99998
PU = solnU$root[2]; PU
## [1] 60
solnE = multiroot(fcnSDE, c(1, 1))
QE = solnE$root[1]; QE
## [1] 90
PE = solnE$root[2]; PE
## [1] 20
solnw = multiroot(fcnESEDw, c(1, 1))
QW = solnw$root[1]; QW
## [1] 75
PW1 = solnw$root[2]; PW1
## [1] 40
# ******** Welfare Analysis ****

ChPSU   = A0U*(PW - PU_Au) + (A1U/2)*(PW**2 - PU_Au**2); ChPSU  # change in producer surplus for country U
## [1] -1750
ChCSU   = B0U*(PU_Au - PW) + (B1U/2)*(PU_Au**2 - PW**2); ChCSU  # change in consumer surplus for country U
## [1] 2500
ChPSE   = A0E*(PW - PE_Au) + (A1E/2)*(PW**2 - PE_Au**2); ChPSE  # change in producer surplus for country E
## [1] 2200
ChCSE   = B0E*(PE_Au - PW) + (B1E/2)*(PE_Au**2 - PW**2); ChCSE  # change in consumer surplus for country E
## [1] -1450
NChU = ChPSU + ChCSU; NChU    # net welfare gains for U
## [1] 750
NChE = ChPSE + ChCSE; NChE    # net welfare gains for E
## [1] 750
NChW  = NChU + NChE; NChW     # net gains for the world
## [1] 1500
rm(list=ls()) # will remove ALL objects


# ******** Analytically solving for autarky and free traded equilibriums ****

# Parameters 

A0U = 25     # supply intercept country U
A1U = 1.25   # supply slope country U
B0U = 250    # demand intercept country U
B1U = -2.5   # demand slope country U

A0E = 50    # supply intercept country E
A1E = 2     # supply slope country E
B0E = 125   # demand intercept country E
B1E = -1.75 # demand slope country E

Tariff = 10


PU_Au = (-A0U + B0U)/(A1U - B1U); PU_Au # Autarky price for country U
## [1] 60
PE_Au = (-A0E + B0E)/(A1E - B1E); PE_Au # Autarky price for country E
## [1] 20
QU_Au = (-A0U*B1U + A1U*B0U)/(A1U - B1U); QU_Au # Autarky Quantity for country U
## [1] 100
QE_Au = (-A0E*B1E + A1E*B0E)/(A1E - B1E); QE_Au # Autarky Quantity for country E
## [1] 90
PW = ((A0U - B0U) + (A0E - B0E))/((-A1U + B1U) + (-A1E + B1E)); PW # World Price
## [1] 40
QE_ES = ((A0E - B0E - B0U + A0U)*(A1E - B1E) + (A0E - B0E)*
           (-A1U + B1U - A1E + B1E))/((-A1U + B1U) + (-A1E + B1E)); QE_ES # Excess supply of country U (Note ES=ED)
## [1] 75
QU_ED =  QE_ES # Excess demand of country E (Note ES=ED)

# ******** Numerically solving for autarky and free traded equilibriums ****

fcnSDU = function(x){
  QU = x[1]
  PU = x[2]
  
  SupplyU = QU - (A0U + A1U*(PW + Tariff)) # QU - (A0U + A1U*PU) = 0 in equilibrium
  DemandU = QU - (B0U + B1U*(PW + Tariff)) # QU - (B0U + B1U*PU) = 0 in equilibrium
  return(c(SupplyU,DemandU))
}

fcnSDE = function(x){
  QE = x[1]
  PE = x[2]
  
  SupplyE = QE - (A0E + A1E*PE) # QE - (A0E + A1E*PE) = 0 in equilibrium
  DemandE = QE - (B0E + B1E*PE) # QE - (B0E + B1E*PE) = 0 in equilibrium
  return(c(SupplyE,DemandE))
}

fcnESEDw = function(x){
  QW = x[1]
  PW = x[2]
  
  ESupplyw = QW - ((B0U + B1U*PW) - (A0U + A1U*PW))
  EDemandw = QW - ((A0E + A1E*PW) - (B0E + B1E*PW))
  
  return(c(ESupplyw,EDemandw))
}

library("rootSolve")
solnU = multiroot(fcnSDU, c(1, 1))
## Warning in stode(y, times, func, parms = parms, ...): error during
## factorisation of matrix (dgefa); singular matrix
## diagonal element is zero 
## [1] 2
## Warning in stode(y, times, func, parms = parms, ...): steady-state not
## reached
QU = solnU$root[1]; QU
## [1] 1
PU = solnU$root[2]; PU
## [1] 1
solnE = multiroot(fcnSDE, c(1, 1))
QE = solnE$root[1]; QE
## [1] 90
PE = solnE$root[2]; PE
## [1] 20
solnw = multiroot(fcnESEDw, c(1, 1))
QW = solnw$root[1]; QW
## [1] 75
PW1 = solnw$root[2]; PW1
## [1] 40
# ******** Welfare Analysis ****

ChPSU   = A0U*(PW - PU_Au) + (A1U/2)*(PW**2 - PU_Au**2); ChPSU  # change in producer surplus for country U
## [1] -1750
ChCSU   = B0U*(PU_Au - PW) + (B1U/2)*(PU_Au**2 - PW**2); ChCSU  # change in consumer surplus for country U
## [1] 2500
ChPSE   = A0E*(PW - PE_Au) + (A1E/2)*(PW**2 - PE_Au**2); ChPSE  # change in producer surplus for country E
## [1] 2200
ChCSE   = B0E*(PE_Au - PW) + (B1E/2)*(PE_Au**2 - PW**2); ChCSE  # change in consumer surplus for country E
## [1] -1450
NChU = ChPSU + ChCSU; NChU    # net welfare gains for U
## [1] 750
NChE = ChPSE + ChCSE; NChE    # net welfare gains for E
## [1] 750
NChW  = NChU + NChE; NChW     # net gains for the world
## [1] 1500