Useful Background: Proposition 9

Existence of RCESE given separability and other restrictions

Suppose every individual \(i\) has a consumption utility function \(v^i(\cdot)\) of the form \[\begin{align} v^i(x^i) = \sum_{\ell=1}^{L-1} \tilde{v}^i_\ell(x^i_\ell) + t^i \cdot x^i, \end{align}\] where the coefficient of relative risk aversion of each \(\tilde{v}^i_\ell\) is everywhere within \((0,1]\). Suppose also that, for all \(i\) and all \(\ell < L\), \(\tilde{v}^i_\ell\) satisfies the lower and upper Inada conditions, and \[\begin{align} t^i_\ell & \geq \bar{\bar{\psi}}^i_\ell, \\ t^i_L & > \bar{P}_\ell\Big(t^i_\ell + \bar{\bar{\psi}}^i_\ell\Big), \text{ and} \\ \omega^i_\ell & > \tilde{v}^{i\prime-1}_\ell\Big(\frac{t^i_L}{\bar{P}_\ell}-t^i_\ell-\bar{\bar{\psi}}^i_\ell\Big). \end{align}\] Then a RCESE exists. Furthermore, in any such RCESE, \(\psi^i_L = 0 \; \forall i\).

Example

Setup

We will work in a two goods case. With two types.Type \(i\) cares about supply and type \(j\) is fully selfish and only cares about personal consumption.

Endowment and other definitions

Let \(N\) be population and set \(\omega^i_\ell=\omega^j_\ell=1,\forall\ell\).

Each type has \(N/2\) individuals. \(N\) is only meaningfully even.

\[ \omega=\frac{N}{2}(\omega^i_1,\omega^i_2)+\frac{N}{2}(\omega^j_1,\omega^j_2)=(N,N) \]

# N <- 1e+06

# N <-1 # Not even so it wont work when there's two types
# N <- 2
# N <- 1e+1
# N <- 1e+2
# N <- 1e+4
N <- 1e+6
# Change to taste!

ag_endowment <- c(N,N)
indiviudal_i_endowment <- c(1,1)
indiviudal_j_endowment <- c(1,1)

message(paste("N has been set to",N))
N has been set to 1e+06

From Proposition 9, \[\psi^i_2 = \psi^j_2= 0.\] We also define from the outset the externality preferences over supply

\[w^i(s)=s_1+2s_2,\] \[w^j(s)=0.\]

psi2_i=0
psi2_j=0


# To automate derivatives we also type it as an expression

w_i_expression = makeFun(s1+2*s2 ~ s1 & s2) #can edit w's expression here, leave '~ s1 & s2' untouched

w_i_s1 = D(w_i_expression(s1,s2) ~ s1) #partial derivative wrt s1 for later
w_i_s1
## function (s1, s2) 
## 1
w_i_s2 =  D(w_i_expression(s1,s2) ~ s2) #partial derivative wrt s2 for later
w_i_s2
## function (s1, s2) 
## 2
# Same for type j
w_j_expression = makeFun(0 ~ s1 & s2) #can edit w's expression here, leave '~ s1 & s2' untouched

w_j_s1 = D(w_j_expression(s1,s2) ~ s1) #partial derivative wrt s1 for later
w_j_s2 =  D(w_j_expression(s1,s2) ~ s2) #partial derivative wrt s2 for later




externalities_function=w_i_expression

Adjusting supply so that it scales with population size

Equation should be:

\[\frac{s_2}{N} = 2 - \frac{1}{2- \dfrac{s_1}{N}}\] \[\Rightarrow s_2 = 2N - \frac{N^2}{2N - s_1} \]

f_ppf =function(x) {2*N-(N^2)/(2*N-x)}

base <-
  ggplot() +
  xlim(-0.3, 1.8*N)+
  ylim(-0.3,2*N)

base + geom_function(fun = f_ppf)+ggtitle(paste("PPF with N =",N)) +
  xlab("s_1") + ylab("s_2")+geom_hline(yintercept = 0)+geom_vline(xintercept = 0)

Finding \(J_s(p)\)

  • Supply function should be a function of the prices, such that suggested (scaled) Production Possibilities Frontier (PPF) becomes: \(s_2 = 2N - \frac{N^2}{2N - s_1}\)

Finding aggregate production \(y(p)\)

Production \(y(p)\) already factors in inputs and outputs. The firm’s ``problem” (considering total output, not firm’s production) is to maximise profit \(p\cdot y(p)\)

\[ \max_{y_1,y_2} \quad p\cdot y(p) \]

Substituting \(y_2\) and using the definition of \(s_i \equiv \omega_i +y_i(p), \forall i\), the problem becomes

\[ \max_{y_1} \quad p_1y_1 + p_2 \left(N -\frac{N^2}{N - y_1}\right) \]

FOC:

\[\begin{align*} &p_1 - p_2 \frac{N^2}{(N - y_1)^2} = 0 \\ &\Rightarrow (N - y_1)^2 = N^2\frac{p_2}{p_1} \end{align*}\]

Hence \[\begin{equation} y_1(p) = N-N\left(\frac{p_2}{p_1}\right)^{1/2}. \label{y1} \end{equation}\]

And by symmetry, \[\begin{equation} \label{y2} y_2(p) = N-N\left(\frac{p_1}{p_2}\right)^{1/2}. \end{equation}\]

[Notice that alternatively \(y_1(p) = \frac{\left(p_{1}+\sqrt{p_{1} p_{2}}\right) N}{p_{1}}, y_2(p) = \frac{\left(p_{2}+\sqrt{p_{1} p_{2}}\right) N}{p_{2}}\)but these solutions do not sit on the PPF we’re interested in.]

Finding aggregate supply \(s(p)\)

\(s(p) \equiv \omega +y(p)\) is hence \[ s(p) = (y_1(p) + N, y_2(p) + N)=N \left(2-\left(\frac{p_2}{p_1}\right)^{1/2}, 2-\left(\frac{p_1}{p_2}\right)^{1/2}\right) \]

supply <- function(p1,p2){
  return(N*c(2-(p2/p1)^(1/2),2-(p1/p2)^(1/2)))
}

Furthermore, due to the non-negativity constraint on aggregate supply, we have \[\begin{equation} \frac{p_2}{4} \le p_1 \le 4p_2 \end{equation}\]

Or in the price ratio form

\[\begin{equation} \frac{1}{4} \le \frac{p_1}{p_2} \le 4 \end{equation}\]

and

\[\begin{equation} \frac{1}{4} \le \frac{p_2}{p_1} \le 4. \end{equation}\] 
# Maximum price ratio (based on notes)
bound_priceratio <- 4

We then have that

\[\begin{equation*} \begin{split} J_s &= \begin{pmatrix} \dfrac{\partial s_1(p)}{\partial p_1} & \dfrac{\partial s_1(p)}{\partial p_2} \\ \dfrac{\partial s_2(p)}{\partial p_1} & \dfrac{\partial s_2(p)}{\partial p_2} \end{pmatrix} \\ &= \begin{pmatrix} \dfrac{N}{2}\left(\dfrac{p_2}{p_1}\right)^{-1/2}\dfrac{p_2}{p_1^2} & -\dfrac{N}{2}\left(\dfrac{p_2}{p_1}\right)^{-1/2}\dfrac{1}{p_1} \\ -\dfrac{N}{2}\left(\dfrac{p_1}{p_2}\right)^{-1/2}\dfrac{1}{p_2} & \dfrac{N}{2}\left(\dfrac{p_1}{p_2}\right)^{-1/2}\dfrac{p_1}{p_2^2} \end{pmatrix} \\ \end{split} \end{equation*}\] Or \[\begin{equation} \begin{split} J_s(p) & = \frac{N}{2} \begin{pmatrix} p_2^\frac{1}{2} p_1^{-\frac{3}{2}} & -(p_1p_2)^{-\frac{1}{2}}\\ -(p_1p_2)^{-\frac{1}{2}} & p_1^\frac{1}{2}p_2^{-\frac{3}{2}} \end{pmatrix} \end{split} \end{equation}\]

J_s <- function(p1,p2){
  matrix(
    N/2*c(p2^(1/2)*p1^(-3/2),
         -(p1*p2)^(-1/2),
         -(p1*p2)^(-1/2),
          p1^(1/2)*p2^(-3/2)
    ),
    nrow = 2 # nrow: number of rows (for a matrix)
  )
}

Ethical Weights’ Bounds

Based on our supply function and externality preferences, we can eventually construct \(\bar{\bar{\psi}}\).

We had \[w^i(s)=-s_1,\] \[w^j(s)=0.\]

Hence we obtain \[\nabla^T w^i(s(p)) =\left( \dfrac{\partial w^i(s)}{\partial s_1} , \dfrac{\partial w^i(s)}{\partial s_2}\right)= \left(\frac{s_1}{((s_1)^2+(s_2)^2)^\frac{1}{2}},\frac{s_2}{((s_1)^2+(s_2)^2)^\frac{1}{2}}\right)\] and \[\nabla^T w^j(s(p)) =\left( \dfrac{\partial w^i(s)}{\partial s_1} , \dfrac{\partial w^i(s)}{\partial s_2}\right)= \left(0,0\right).\]

For the time being let us focus on type \(i\).

Bounds for type \(i\)

# Gradient of w:


# [Exact formulae commented out]
# Gradient of w:
grad_externalities <- function (p1,p2){
  s1 = supply(p1,p2)[1]
  s2 = supply(p1,p2)[2]
  # den = (s1^2+s2^2)^(1/2)
  # 
  # return(c(s1/den,s2/den))
  return(c(w_i_s1(s1,s2),w_i_s2(s1,s2)))
}

# matrix version
grad_externalities2 <- function (p1,p2){
  s1 = supply(p1,p2)[1]
  s2 = supply(p1,p2)[2]
  # den = (s1^2+s2^2)^(1/2)
  # 
  # return(matrix(c(s1/den,s2/den), nrow = 1))
  return(matrix(c(w_i_s1(s1,s2),w_i_s2(s1,s2)), nrow = 1))
}


# some tests
# grad_externalities(1,0.25000000000001)
# grad_externalities2(1,0.25000000000001)
# grad_externalities(1,0.2501)
# grad_externalities2(1,0.2501)
grad_externalities(1,0.25)
## [1] 1 2
grad_externalities2(1,0.25)
##      [,1] [,2]
## [1,]    1    2
grad_externalities(1,4)
## [1] 1 2
grad_externalities2(1,4)
##      [,1] [,2]
## [1,]    1    2

Prop 9 imposes that the coefficient of relative risk aversion is in \((0,1]\). Taking the case where \(\gamma\) equals 1, let \[\begin{equation} \label{logutil} \tilde{v}^i_1 =\tilde{v}^j_1 \equiv \ln(x_1) \end{equation}\]

The code can be used to find \(t^i_1\) and \(t^i_2\) for any given price level using the fact that:

\[\begin{align} \bar{\bar{\psi}}^i_\ell &\equiv \max_{p \in P, k\leq L, m\leq L} \big|J_s(p) e_k \cdot \nabla w^i(s(p))\big|\big|e_\ell \cdot J_s(p) e_m\big|^{-1} \;\;\; (\ell < L).\\ t^i_\ell &\geq \bar{\bar{\psi}}^i_\ell \;\;\; \forall i, \; \forall \ell < L \label{maxrescale} \end{align}\] 
bound_psi1 <- function(p1,p2){
# Take a good l, we want to construct barbarpsi_l
# we only need to do this for 1 here because it's the only l\leq L=2
e_l=c(1,0)
# create a matrix of m times k (2X2) st each row is the kth value and m corresponds to col
boundcandidates=matrix(numeric(4), nrow=2)

e_m=c(0,0)
for (k in 1:2){
  # vary the position of the 1 in e_k
  e_k=c(0,0)
  e_k[k]=1
    # numerator= absolute value of  e_k'  J_s(p)' gradient of externalities at p (for varying k \leq L=2)
  num= abs( t(e_k) %*% t(J_s(p1, p2)) %*% grad_externalities(p1,p2))
  
  for (m in 1:2) {
    # vary the position of the 1 in e_m
    e_m=c(0,0)
    e_m[m]=1
    # denominator= absolute value of e_l' J_s(p) e_m  (for varying m \leq L=2)
    denom= abs(t(e_l)  %*% J_s(p1, p2) %*% e_m)
# assign to matrix
    boundcandidates[k,m]=num/denom
  }
}

# we can then get the max candidate:
return(max(boundcandidates))
}

We then need to loop this through possible prices as well. As we do this we must construct ts such that Proposition 9’s conditions hold. Namely,

\[\begin{align} t^i_1 & \geq \bar{\bar{\psi}}^i_1, \\ t^i_2 & > \bar{P}_1\Big(t^i_1 + \bar{\bar{\psi}}^i_1\Big), \\ \end{align}\] and \[\begin{align} \omega^i_1 &> \tilde{v}^{i\prime-1}_1\Big(\frac{t^i_2}{\bar{P}_1} - t^i_1 - \bar{\psi}^i_1\Big)\\ \iff t^i_2 &> \bar{P}_1 \big( \tilde{v}^{i\prime}_1(\omega^i_1) + t^i_1 + \bar{\psi}^i_1 \big). \end{align}\]

Where the inverse of the first derivative of \(\ln(\cdot)\) is

\[\tilde{v}^{i\prime-1}(x)=\frac{1}{x}\]

and where we recall that \(\omega^i_1=1\) such that \(\tilde{v}^{i\prime}_1(\omega^i_1)=1\).

# ts in selfish utility function
t_s <- function (p1,p2) {
  c(bound_psi1(p1,p2) + 1, bound_priceratio*(1+2*bound_psi1(p1,p2))+1)
# seems like it could also be c(bound_psi1(p1,p2), bound_priceratio*(1+2*bound_psi1(p1,p2))+1)
}

t1=0
t2=0

p1=1

for (p2 in seq(from = 0.25, to = 4, by = 0.05)) {
  # print(t_s(p1,p2))
  
  if (t_s(p1,p2)[1]>t1) {
    t1=ceiling(t_s(p1,p2)[1])
  }
  if (t_s(p1,p2)[2]>t2) {
    t2=ceiling(t_s(p1,p2)[2])
  }
  #cycle through prices to find bound of psis and therefore ts
}
# to set manually uncomment these two lines below
t1=1/3
t2=1
# t2=67/90 + sqrt(73)/18
# t1=1/3 #case 1/3,1,0,1 gives two equilibria
# t2=1 #case 1/3,1,0,1 gives two equilibria

Bounds for type \(j\)

We can again use

\[\begin{align} \bar{\bar{\psi}}^j_\ell &\equiv \max_{p \in P, k\leq L, m\leq L} \big|J_s(p) e_k \cdot \nabla w^j(s(p))\big|\big|e_\ell \cdot J_s(p) e_m\big|^{-1} \;\;\; (\ell < L).\\ t^j_\ell &\geq \bar{\bar{\psi}}^j_\ell \;\;\; \forall i, \; \forall \ell < L. \end{align}\]

In fact from the paper we defined \[\psi^j(p,x(\cdot)) \equiv \big(J_s(p)\, F(p,x(\cdot))\big)^T \nabla w^j(s(p)).\]

Here, \[\nabla^T w^j(s(p)) = \left(0,0\right).\]

So it immediately follows that

\[ \bar{\bar{\psi}}^j_\ell= \psi^j_\ell=0, \;\;\; \forall \ell.\] We construct \(t\)s such that Proposition 9’s conditions hold. Namely,

\[\begin{align} t^j_1 & \geq \bar{\bar{\psi}}^j_1=0, \\ t^j_2 & > \bar{P}_1\Big(t^j_1 + \bar{\bar{\psi}}^j_1\Big)= \bar{P}_1(t^j_1), \\ \end{align}\] and \[\begin{align} \omega^j_1 &> \tilde{v}^{j\prime-1}_1\Big(\frac{t^j_2}{\bar{P}_1} - t^j_1 - \bar{\psi}^j_1\Big)\\ \iff t^j_2 &> \bar{P}_1 \big( \tilde{v}^{j\prime}_1(\omega^j_1) + t^j_1 + \bar{\psi}^j_1 \big). \end{align}\]

Where the inverse of the first derivative of \(\ln(\cdot)\) is

\[\tilde{v}^{j\prime-1}(x)=\frac{1}{x}\]

and where we recall that \(\omega^j_1=1\) such that \(\tilde{v}^{j\prime}_1(\omega^j_1)=1\).

In particular, \(t^j_1=0\) and \(t^j_2 = 5 >{\bar{P}_1}\) fulfil the above conditions.

\[\begin{align} t^j_1=0 & \geq \bar{\bar{\psi}}^j_1=0, \\ t^j_2 =5 & > \bar{P}_1\Big(t^j_1 + \bar{\bar{\psi}}^j_1\Big)= \bar{P}_1(t^j_1)=4(0), \\ \end{align}\] and \[\begin{align} t^j_2=5 & > \bar{P}_1 \big( \tilde{v}^{j\prime}_1(\omega^j_1) + t^j_1 + \bar{\psi}^j_1 \big)=4(1+0+0). \end{align}\]

# t1_j=0
# t2_j=5
# t1_j=t1
# t2_j=t2
t1_j=0 #case 1/3,1,0,1 gives two equilibria
t2_j=t2#case 1/3,1,0,1 gives two equilibria
# cutting slack as much as possible t_j would be
# t1_j=0
# t2_j=bound_priceratio+1
  c(t1_j,t2_j)
## [1] 0 1
psi1_j=0
# psi2_j=0 already above
psi2=psi2_i # now that j is known to give zeros drop the i for psi2_i easier notation

Create implicit demand functions

Individual problem is pinned down by Prop 4 \[\begin{equation} \bar{x}^i(\bar{p},\bar{\pi}) = \max\limits_{x^i} \Big(v^i(x^i) + \psi^i(\bar{p},\bar{x}(\cdot)) \cdot x^i \Big) \;\; \Big| \;\; \bar{p} \cdot x^i \leq \bar{p} \cdot \omega^i + \theta^i \bar{\pi} \;\; \forall i. \end{equation}\] And recalling the form of the utility function given in Prop 9 \[\begin{align} v^i(x^i) = \sum_{\ell=1}^{L-1} \tilde{v}^i_\ell(x^i_\ell) + t^i \cdot x^i= \tilde{v}^i_1(x^i_1) + t_1^ix_1^i+t_2^ix_2^i=\ln(x_1)+ t_1^ix_1^i+t_2^ix_2^i. \end{align}\] with the last two equalities holding in the two goods case with the log utility function. An individual therefore maximises: \[ \max_{x_1^i, x_2^i} \quad \ln(x_1^i) + (t_1^i+ \psi_1^i) x_1^i + (t_2^i + \psi_2^i) x_2^i \] \[ \text{ s.t. } p_1x_1^i + p_2x_2^i = p_1 \omega_1^i + p_2 \omega_2^i + \frac{\pi(p)}{N} \] We have that $t^i_1=

FOC (assuming that both goods will be consumed in positive quantities): \[ \frac{1}{p_1}\left(\frac{1}{x_1^i} + t_1^i + \psi_1^i\right) = \frac{t_2^i + \psi_2^i}{p_2} \] \[ \Rightarrow \frac{1}{x_1^i} = \frac{p_1}{p_2}(t_2^i+\psi_2^i) - (t_1^i + \psi_1^i) \] \[\begin{equation} \label{x_1} \chi_1^i(p) = \frac{p_2}{p_1(t_2^i+\psi_2^i) -p_2(t_1^i + \psi_1^i)} \end{equation}\]

Imposing budget, \[ \chi_2^i(p) = \frac{p_1}{p_2}\left(\omega_1^i - \frac{p_2}{p_1(t_2^i+\psi_2^i) - p_2(t_1^i + \psi_1^i)}\right) + \omega_2^i + \frac{\pi(p)}{Np_2} \]

We can use \(\pi= p\cdot y(p)=p_1y_1+p_2y_2\). By equations on output, \[\pi= p_1(N-N\left(\frac{p_2}{p_1}\right)^{1/2})+ p_2(N-N\left(\frac{p_1}{p_2}\right)^{1/2}),\] such that \[\begin{align*} \frac{\pi}{Np_2}=& \frac{p_1}{p_2}(1-\left(\frac{p_2}{p_1}\right)^{1/2})+ (1-\left(\frac{p_1}{p_2}\right)^{1/2})\\ =&\frac{p_1}{p_2}-2\left(\frac{p_1}{p_2}\right)^{1/2}+1\\ =&\left(\left(\frac{p_1}{p_2}\right)^{1/2}-1\right)^2. \end{align*}\]

Which we substitute back into \(\chi_2^i(p)\),

\[\begin{equation} \label{x_2} \chi_2^i(p) = \frac{p_1}{p_2}\left(\omega_1^i - \frac{p_2}{p_1(t_2^i+\psi_2^i) - p_2(t_1^i + \psi_1^i)}\right) + \omega_2^i + \left(\left(\frac{p_1}{p_2}\right)^\frac{1}{2}-1\right)^2. \end{equation}\] 
demand <- function(p1,p2,t1,t2,psi1,psi2) {
  #Some intermediate terms------
  r=p1/p2
  a=(t2+psi2)
  b=(t1+psi1)
  d=p1*a-p2*b
  #------
  return((N/2)*c(p2/d,r*(1-p2/d)+1+(r^(1/2)-1)^2))
  # used to be for only one individual and not have N  
  #N/2 since now it's for half of the population, i.e. aggregate for each type
}

excessdemand <- function(supply,demand){
  return(demand-supply)
}

Notice that we did not use anything specific about \(i\) above, hence all the results in this section hold for \(j\) in the same way, just with different \(t\)s and \(\psi\)s. The symmetry between \(i\) and \(j\) carries in the sections below.

Construct the matrix \(J_z(p)\)

By definition, \(J_z = J_{\chi} - J_s\)

Taking the derivative with respect to prices of the expressions above:

\[\begin{align} \frac{\partial \chi_1^i(p)}{\partial p_1} = -\frac{p_2(t_2^i+\psi_2^i)}{[p_1(t_2^i+\psi_2^i) -p_2(t_1^i + \psi_1^i)]^2} \end{align}\] \[\begin{align} \frac{\partial \chi_1^i(p)}{\partial p_2} = \frac{1}{p_1(t_2^i+\psi_2^i) -p_2(t_1^i + \psi_1^i)} + \frac{p_2(t_1^i +\psi_1^i)}{[p_1(t_2^i+\psi_2^i) -p_2(t_1^i + \psi_1^i)]^2} \end{align}\] \[\begin{equation} \begin{aligned} \frac{\partial \chi_2^i(p)}{\partial p_1} =& \frac{1}{p_2}\omega_1^i - \frac{1}{p_1(t_2^i+\psi_2^i) - p_2(t_1^i + \psi_1^i)} + \frac{p_1(t_2^i+\psi_2^i)}{[p_1(t_2^i+\psi_2^i) - p_2(t_1^i + \psi_1^i)]^2} \\ &+\frac{1}{p_2}\left[1-\left(\frac{p_2}{p_1}\right)^{1/2}\right] \end{aligned} \end{equation}\] \[\begin{equation} \frac{\partial \chi_2^i}{\partial p_2} = -\frac{p_1}{p_2^2}\omega_1^i - \frac{p_1(t_1^i+\psi_1^i)}{[p_1(t_2^i+\psi_2^i) - p_2(t_1^i + \psi_1^i)]^2} - \frac{1}{p_2}\left[\frac{p_1}{p_2} -\left(\frac{p_1}{p_2}\right)^{1/2}\right] \end{equation}\]

Putting it all together, the Jacobian taking into account all members of type \(k\),

\[\begin{equation} \begin{split} J_{\chi^k}(p) &= \frac{N}{2} \begin{pmatrix} \dfrac{\partial \chi_1^k(p)}{\partial p_1} & \dfrac{\partial \chi_1^k(p)}{\partial p_2} \\ \dfrac{\partial \chi_2^k(p)}{\partial p_1} & \dfrac{\partial \chi_2^k(p)}{\partial p_2} \end{pmatrix}, \end{split} \end{equation}\]

with each entry matching the 4 equations above (where \(i\) is in fact any type \(k\)).

##Jacobian of Implicit Demand------

J_x <- function(p1,p2,t1,t2,psi1,psi2) {
  Df <- matrix(numeric(4),nrow = 2,ncol = 2) #creates matrix of 4 zeros organised as entries into nrow rows and ncol cols
  #Some intermediate terms------
  r=p1/p2
  a=(t2+psi2)
  b=(t1+psi1)
  d=p1*a-p2*b
  #------
  Df[1,1] <- -p2*a/(d^2)
  Df[1,2] <- 1/d+(p2*b)/(d^2)
  Df[2,1] <- 1/p2-1/d+p1*a/(d^2)+(1/p2)*(1-r^(-1/2))
  Df[2,2] <- -p1/(p2^2)-(p1*b)/(d^2)-(1/p2)*(r-r^(1/2))
  return((N/2)*Df)
}
J_x(1,1,0,1,0,0)
##        [,1]   [,2]
## [1,] -5e+05  5e+05
## [2,]  5e+05 -5e+05

(Testing analytic example)

J_x(1,1,0,1,0,0)#J_x for j which is what i best responds to when calculating -Jz leading to F
##        [,1]   [,2]
## [1,] -5e+05  5e+05
## [2,]  5e+05 -5e+05
J_s(1,1)
##        [,1]   [,2]
## [1,]  5e+05 -5e+05
## [2,] -5e+05  5e+05
c(t1,t2,psi2)
## [1] 0.3333333 1.0000000 0.0000000
J_x(1,1,t1,t2,0,0)
##          [,1]    [,2]
## [1,] -1125000 1125000
## [2,]   875000 -875000
Ftest=matrix(c(1/3,0,0,0),nrow=2)
Ftest
##           [,1] [,2]
## [1,] 0.3333333    0
## [2,] 0.0000000    0
# AA^gA=A for A^g the generalised inverse of A
Jac_test=J_s(1,1)-J_x(1,1,0,1,0,0)
Jac_test
##        [,1]   [,2]
## [1,]  1e+06 -1e+06
## [2,] -1e+06  1e+06
Jac_test %*% Ftest %*% Jac_test
##               [,1]          [,2]
## [1,]  333333333333 -333333333333
## [2,] -333333333333  333333333333

Moreover, recall \[\begin{equation} \begin{split} J_s(p) & = \frac{N}{2} \begin{pmatrix} p_2^\frac{1}{2} p_1^{-\frac{3}{2}} & -(p_1p_2)^{-\frac{1}{2}}\\ -(p_1p_2)^{-\frac{1}{2}} & p_1^\frac{1}{2}p_2^{-\frac{3}{2}} \end{pmatrix}. \end{split} \end{equation}\]

We can hence write \(J_z = J_{\chi} - J_s\).

Obtaining \(F\)

\[\begin{align} F(p,x(\cdot)) \; \triangleq \; &\text{The generalized inverse $F$ of $-J_z(p)$, where $z(\cdot)$ is implied by} \nonumber\\ &\text{$(\mathcal{E}, x(\cdot))$, with $F\delta(p,x(\cdot)) = \mathbb{0}_L$ and whose bottom row equals }\mathbb{0}_L^T \nonumber \\ \end{align}\]

And here:

\(\bar{\delta} = e_2/\bar{p}_2\) where \(e_2=\begin{pmatrix} 0\\ 1 \end{pmatrix}\)

\[\begin{align} \psi^i(p,x(\cdot)) \; \triangleq \; & \big(J_s(p)\, F(p,x(\cdot))\big)^T \nabla w^i(s(p)) \;\;\; \\ \psi(p,x(\cdot)) \; \triangleq \; & \text{The $L \times I$ matrix with column $i$ equal to $\psi^i(p,x(\cdot))$} \nonumber \end{align}\]

Singular Value Decomposition

We follow the construction in the paper Appendix A4.

We will now construct the unique generalised inverse \(F\) of \(-J_z(\bar{p})\) such that

Let \(U N V^T\) denote a singular value decomposition of \(J_z(\bar{p})\) with \(N_{1 1} \neq 0\) and \(N_{2 2} = 0\). Because \(F\) is a generalised inverse of \(J_z(\bar{p})\), we must by (107) have \[\begin{align} F = V\bigg[ \begin{matrix} N_1^{-1} & A \\ B & C \end{matrix} \bigg]U^T, \end{align}\] where \(N_1\) is the \((1) \times (1)\) principal submatrix of \(N\) (ie scalar) and \(A\), \(B\), and \(C\) are all scalar too.

Since \(V\) is invertible, (109) reduces to \[\begin{align} \bigg[ \begin{matrix} N_1^{-1} & A \\ B & C \end{matrix} \bigg] \bar{\bar{\delta}} &= \begin{pmatrix} 0\\ 0 \end{pmatrix} \\ \implies A &= -\frac{\bar{\bar{\delta}}_1}{\bar{\bar{\delta}}_2} \, \frac{1}{N_{1 1}}, \;\; \text{ and} \\ C &= -\frac{\bar{\bar{\delta}}_1}{\bar{\bar{\delta}}_2} B , \label{C} \end{align}\] where \(\bar{\bar{\delta}} \triangleq U^T \bar{\delta}\).

And here:

\(\bar{\delta} = e_2/\bar{p}_2\) where \(e_2=\begin{pmatrix} 0\\ 1 \end{pmatrix}\), giving

\[\bar{\delta} = \begin{pmatrix} 0\\ \frac{1}{\bar{p}_2} \end{pmatrix}.\]

We will now impose the constraint that the bottom row of \(F\) consists of zeroes. Since \(U^T\) is invertible, when the bottom row of \(F\) consists of zeroes,

\[\begin{align} B &= - \frac{V_{21}}{V_{22}} \, \frac{1}{N_{1 1}}. \end{align}\] This also gives us \(C\). \(F\) is thus fully constructed. Following the notation in the code,

\[\begin{align} F = V\bigg[ \begin{matrix} N_1^{-1} & A \\ B & C \end{matrix} \bigg]U^T=V N_{SVD} U^T. \end{align}\] 
F <- function(J_s,J_x) {
  X=-(J_x-J_s)
  # X=scale(X)
  (s <- svd(X)) #crucial step, SVD happening here!
  U=s$u
  s$d[length(s$d)]=0 #set last one to zero, should be close to zero already anyway
  V=s$v
  D <- diag(s$d) #what the paper calls N
  # s$u %*% D %*% t(s$v) # X = U D V'
  N_1_inv= 1/D[1,1] #what the paper calls N_1^{-1}. note solve(D) gives inverse of D
  
  delta_bar=matrix(c(0,1/p2))
  # for now we imagine p2 is really the eqm price (because all we need is for it to hold true in the eqm prices we're trying to find)
  
  delta_barbar= t(U) %*% delta_bar
  
  A=-(delta_barbar[1]/delta_barbar[2])*(1/D[1,1])
  B=-(V[2,1]/V[2,2])*(1/D[1,1])
  C=-(delta_barbar[1]/delta_barbar[2])*B
  # follow equations in notes for the above

  N_SVD=matrix(c(0,0,0,0), nrow = 2)
  N_SVD[1,1]=N_1_inv
  N_SVD[1,2]=A
  N_SVD[2,1]=B
  N_SVD[2,2]=C
  
  return(V %*% N_SVD %*% t(U))
  }

All that is left is to find the fixed point of \(\psi^i(p,x(\cdot))\), or the root of \(\psi^i(p,x(\cdot))-\psi\), which can be done solving the linear system in R. This will give psis for different prices and then, to find the equilibrium we just ensure that excess demand is zero.

Once \(\psi\) is determined, the equilibrium prices that clear the market can be found. Since relative prices are the key, we can set \(p_1=1\) and vary \(p_2\) between \(0.25p_1\) and \(4p_1\).

c(p1,p2)
## [1] 1 4
c(t1,t2,psi1_j,psi2_j)
## [1] 0.3333333 1.0000000 0.0000000 0.0000000
c(t1_j,t2_j,psi1_j,psi2_j)
## [1] 0 1 0 0
fnpsi_minuspsi <- function(psi1) {
  totalJx=J_x(p1,p2,t1,t2,psi1,psi2_i)+J_x(p1,p2,t1_j,t2_j,psi1_j,psi2_j)#we need to use Jx so we add J_x of i and J_x of j
  F=F(J_s(p1,p2),totalJx)
  Js=J_s(p1,p2)
  y <-t(Js %*% F) %*% t(grad_externalities2(p1,p2))
  y[1] <- y[1]-psi1
  # y[2] <- y[2]-psi2
  return(y[1])
}
# original general function (which should apply to when there's only one type)
# fnpsi_minuspsi <- function(psi1) {
#   F=F(J_s(p1,p2),J_x(p1,p2,t1,t2,psi1,psi2))
#   J=J_s(p1,p2)
#   y <-t(F%*%J)%*%t(grad_externalities2(p1,p2))
#   y[1] <- y[1]-psi1
#   # y[2] <- y[2]-psi2
#   return(y[1])
# }
# should only be in terms of psis at this point (if t stuff ran well)
# the above J_x is the Jacobian of j's implicit demand to which i is best responding


# nleqslv(psi1start, fnpsi_minuspsi, control=list(btol=.01))



# MARKET CLEARING -------
# does market clear, is z(p)=0?
# if it has tiny almost zero entries then fine! record these prices and then we want to create dataframe with the prices, the price ratio, the demand, the psis, the ts

# function for market clearing

closetozero <- function(z1,z2){
  a=0.001*N
  # a=0.001*N originally 
  # tolerance can be made to depend on size of economy N
  if (z1 < a & z1 > -a & z2 < a & z2 > -a) {
    return("All clear!")
  }
  return(" ")
}

# function for market low loss [comes later]

lowloss <- function(loss){
  # a=300*N
  a=0.0001*N #originally
  # tolerance can be made to depend on size of economy N
  if (loss < a) {
    return("Low loss!")
  }
  return(" ")
}

# --------
# guess value for psi1
psi1start=0.1

# Setup for for loop----
start=0.25
end=4
# step=0.0001 
step=0.0025 # original

 # step=0.0000001 # can change step to tiny a final time
num_of_steps=(end-start)/step+1

table_colnames=c("p1","p2","t1_i","t2_i","psi1_i","psi2_i=psi2_j=psi1_j","t1_j","t2_j","z1","z2","doubleclearing")
#doubleclearing is when Both markets clear +-0.001*N

# Big table with key parameters----
table = matrix(numeric(num_of_steps*length(table_colnames)),nrow = num_of_steps ,ncol = length(table_colnames) )
colnames(table)=table_colnames
# col names of matrix in order

# Now use a For loop to produce the big table elements--------

# Really precise estimates from tests:------------------------
# With step=0.0000001, range [3.35,3.36]
# N p2_equilibrium  loss
# 1 3.3551559000000001909   1.588051845295347325e-18
# 1e+06 3.3528541999999998957   1.4720168453199354205e-06

# previously 3.353535353...=332/99
# ------------------------------------------------------------

for (p2 in seq(from = start, to = end, by = step)) {
  if(p2==start){
    count=1
  }
  psi1=nleqslv(psi1start, fnpsi_minuspsi, control=list(btol=.01))$x #this is the key step finding fixed points!
z=  excessdemand(supply(p1,p2),demand(p1,p2,t1,t2,psi1,psi2)+demand(p1,p2,t1_j,t2_j,psi1_j,psi2_j)) #for total agg dem we add (agg) demand of both types
  
  table[count,1]= p1 #should always be one
  table[count,2]= p2
  table[count,3]= t1 #the ones without subindex are i
  table[count,4]= t2 #the ones without subindex are i
  table[count,5]= psi1
  table[count,6]= psi2 #should always be zero
  
  table[count,7]= t1_j 
  table[count,8]= t2_j

  table[count,9]=  z[1]
  table[count,10]= z[2]
  table[count,11]= closetozero(z[1], z[2])
  
  count=count+1
}

Results

We obtain the table summarising the results.

head(table)
##      p1  p2       t1_i                t2_i psi1_i             
## [1,] "1" "0.25"   "0.333333333333333" "1"  "-4.3053644064728" 
## [2,] "1" "0.2525" "0.333333333333333" "1"  "-4.24688404954674"
## [3,] "1" "0.255"  "0.333333333333333" "1"  "-4.18964229271134"
## [4,] "1" "0.2575" "0.333333333333333" "1"  "-4.13360163746281"
## [5,] "1" "0.26"   "0.333333333333333" "1"  "-4.07872606270429"
## [6,] "1" "0.2625" "0.333333333333333" "1"  "-4.0249809532087" 
##      psi2_i=psi2_j=psi1_j t1_j t2_j z1                  z2                
## [1,] "0"                  "0"  "1"  "-1312280.72615213" "5249122.9046085" 
## [2,] "0"                  "0"  "1"  "-1307755.66255137" "5179230.3467381" 
## [3,] "0"                  "0"  "1"  "-1303239.86373796" "5110744.56367827"
## [4,] "0"                  "0"  "1"  "-1298733.14261171" "5043623.85480275"
## [5,] "0"                  "0"  "1"  "-1294235.31632054" "4977828.13969437"
## [6,] "0"                  "0"  "1"  "-1289746.20611501" "4913318.88043814"
##      doubleclearing
## [1,] " "           
## [2,] " "           
## [3,] " "           
## [4,] " "           
## [5,] " "           
## [6,] " "

As numeric data frame and condensed:

frame <- as.data.frame(table)
frame[, 1:10] <- sapply(frame[, 1:10], as.numeric) #ensure everything but doubleclearing is numeric, make sure to change range if columns are added

print(frame[1:4,], digits = 5)
##   p1     p2    t1_i t2_i  psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j       z1
## 1  1 0.2500 0.33333    1 -4.3054                    0    0    1 -1312281
## 2  1 0.2525 0.33333    1 -4.2469                    0    0    1 -1307756
## 3  1 0.2550 0.33333    1 -4.1896                    0    0    1 -1303240
## 4  1 0.2575 0.33333    1 -4.1336                    0    0    1 -1298733
##        z2 doubleclearing
## 1 5249123               
## 2 5179230               
## 3 5110745               
## 4 5043624

And we can see when markets clear \(\pm 0.001N\) for the equilibria candidates below.

equilibria <- frame[frame$doubleclearing == "All clear!", ]  
print(equilibria, digits = 5)
##     p1 p2    t1_i t2_i   psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j         z1
## 301  1  1 0.33333    1 -0.33333                    0    0    1 8.9965e-06
##              z2 doubleclearing
## 301 -8.9965e-06     All clear!

Other checks

We had \[\begin{align} v^i(x^i) = \tilde{v}^i_1(x^i_1) + t^i_1x^i_1+t^i_2x^i_2 \text{ for } L=2, \end{align}\]

Where

\[\begin{align} \tilde{v}^i_1(x^i_1) = \ln (x^i_1) \end{align}\]

Proposition 9 conditions

The main conditions were

\[\begin{align} t^i_1 & \geq \bar{\bar{\psi}}^i_1, \\ t^i_2 & > \bar{P}_1\Big(t^i_1 + \bar{\bar{\psi}}^i_1\Big), \text{ and} \\ \omega^i_1 & > \tilde{v}^{i\prime-1}_1\Big(\frac{t^i_2}{\bar{P}_1}-t^i_1-\bar{\bar{\psi}}^i_1\Big). \end{align}\]

We check the first, that \(t^i_1 \geq \bar{\bar{\psi}}^i_1\) :

# by inspection
range(frame$psi1)
## [1] -4.30536441  0.04967264
range(frame$t1_i)
## [1] 0.3333333 0.3333333
# test
max(frame$psi1)<max(frame$t1_i)
## [1] TRUE

We check the second, that \(t^i_2 > \bar{P}_1\Big(t^i_1 + \bar{\bar{\psi}}^i_1\Big)\) :

# by inspection
range(frame$t2_i)
## [1] 1 1
bound_priceratio*(max(frame$t1_i)+max(frame$psi1))
## [1] 1.532024
# test
max(frame$t2_i)>bound_priceratio*(max(frame$t1_i)+max(frame$psi1))
## [1] FALSE

Lastly, we supposed that endowments \(\omega^i\) satisfy \[\begin{align} t^i_2 &> \bar{P}_1 \big( \tilde{v}^{i\prime}_1(\omega^i_1) + t^i_1 + \bar{\psi}^i_1 \big)\\ \iff \omega^i_1 &> \tilde{v}^{i\prime-1}_1\Big(\frac{t^i_2}{\bar{P}_1} - t^i_1 - \bar{\psi}^i_1\Big), \label{endowmentcond} \end{align}\]

And we know \(\omega^i_1 =1\).

# log(x) derives into 1/x and inverse of 1/x is itself
# inspect
max(frame$t2_i)
## [1] 1
bound_priceratio*(1/1+max(frame$t1_i)+max(frame$psi1))
## [1] 5.532024
# test
max(frame$t2_i)>bound_priceratio*(1/1+max(frame$t1_i)+max(frame$psi1))
## [1] FALSE

In the inverse form:

# RHS must be less than 1
1/(max(frame$t2_i)/bound_priceratio-max(frame$t1_i)-max(frame$psi1))
## [1] -7.51846
1>1/(max(frame$t2_i)/bound_priceratio-max(frame$t1_i)-max(frame$psi1))
## [1] TRUE

Complete economy at equilibria

Here I build an even bigger data frame with the rest of the economy spelled out.

extra_colnames=c("x1_i","x2_i","s1","s2","y1","y2","v_tilde_i","v_i","w_i","u_i","expend_i-budget_i","x1_j","x2_j","v_tilde_j","v_j","w_j","u_j","expend_j-budget_j")

# useful to explicitly define output
output <- function(p1,p2){
   return(c(supply(p1,p2)[1]-N,supply(p1,p2)[2]-N))
}

# and other utility functions

v_tilde_function <- function(x1){
   return(log(x1))
}

v_function <- function(x1,x2,t1,t2){
   return(v_tilde_function(x1)+x1*t1+x2*t2)
}

# u_function is just v_function+ externalities_function


# Big table with everything----
mega_frame=  matrix(numeric(num_of_steps*length(extra_colnames)),nrow = num_of_steps ,ncol = length(extra_colnames) )
colnames(mega_frame)=extra_colnames


for (p2 in seq(from = start, to = end, by = step)) {
  if(p2==start){
    count=1
  }
  # i:
  x1_i=demand(p1=p1,p2=p2,
              t1=frame$t1_i[count],
              t2=frame$t2_i[count],
              psi1=frame$psi1_i[count],
              psi2=frame$psi2_i[count])[1]
  x2_i=demand(p1=p1,p2=p2,
              t1=frame$t1_i[count],
              t2=frame$t2_i[count],
              psi1=frame$psi1_i[count],
              psi2=frame$psi2_i[count])[2]
  w_i=externalities_function(supply(p1,p2)[1],supply(p1,p2)[2])
  v_i=v_function(x1_i,x2_i,
                 t1=frame$t1_i[count],
                 t2=frame$t2_i[count])
  # j:
  x1_j=demand(p1=p1,p2=p2,
              t1=frame$t1_j[count],
              t2=frame$t2_j[count],
              psi1=0,psi2=0)[1]
  x2_j=demand(p1=p1,p2=p2,
              t1=frame$t1_j[count],
              t2=frame$t2_j[count],
              psi1=0,psi2=0)[2]
  w_j=0
  v_j=v_function(x1_j,x2_j,
                 t1=frame$t1_j[count],
                 t2=frame$t2_j[count])
  # add to table:
  mega_frame[count,1]= x1_i
  mega_frame[count,2]= x2_i
  mega_frame[count,3]= supply(p1,p2)[1]
  mega_frame[count,4]= supply(p1,p2)[2]
  mega_frame[count,5]= output(p1,p2)[1]
  mega_frame[count,6]= output(p1,p2)[2]
  mega_frame[count,7]= v_tilde_function(x1_i)
  mega_frame[count,8]= v_i
  mega_frame[count,9]= w_i
  mega_frame[count,10]= v_i+w_i
  mega_frame[count,11]= x1_i*p1+x2_i*p2-(p1*ag_endowment[1]/2+p2*ag_endowment[2]/2+
    (p1*output(p1,p2)[1]+p2*output(p1,p2)[2])/2) #price of cons=expendi=value of endownment+profit (i.e. this should be zero when things work correctly)

  
  mega_frame[count,12]= x1_j
  mega_frame[count,13]= x2_j
  mega_frame[count,14]= v_tilde_function(x1_j)
  mega_frame[count,15]= v_j
  mega_frame[count,16]= w_j
  mega_frame[count,17]= v_j+w_j
  mega_frame[count,18]= x1_j*p1+x2_j*p2-(p1*ag_endowment[1]/2+p2*ag_endowment[2]/2+
    (p1*output(p1,p2)[1]+p2*output(p1,p2)[2])/2) #price of cons=expendi=value of endownment+profit (i.e. this should be zero when things work correctly)
  
    count=count+1
}



mega_frame <- as.data.frame(mega_frame)
mega_frame <- cbind(frame, mega_frame)
head(mega_frame)
##   p1     p2      t1_i t2_i    psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j       z1
## 1  1 0.2500 0.3333333    1 -4.305364                    0    0    1 -1312281
## 2  1 0.2525 0.3333333    1 -4.246884                    0    0    1 -1307756
## 3  1 0.2550 0.3333333    1 -4.189642                    0    0    1 -1303240
## 4  1 0.2575 0.3333333    1 -4.133602                    0    0    1 -1298733
## 5  1 0.2600 0.3333333    1 -4.078726                    0    0    1 -1294235
## 6  1 0.2625 0.3333333    1 -4.024981                    0    0    1 -1289746
##        z2 doubleclearing     x1_i    x2_i      s1       s2       y1         y2
## 1 5249123                62719.27 2749123 1500000     0.00 500000.0 -1000000.0
## 2 5179230                63500.56 2718834 1497506  9925.62 497506.2  -990074.4
## 3 5110745                64284.89 2689176 1495025 19704.91 495024.8  -980295.1
## 4 5043624                65072.28 2660129 1492555 29341.44 492555.4  -970658.6
## 5 4977828                65862.73 2631674 1490098 38838.65 490098.0  -961161.4
## 6 4913319                66656.26 2603795 1487652 48199.85 487652.5  -951800.1
##   v_tilde_i     v_i     w_i     u_i expend_i-budget_i   x1_j    x2_j v_tilde_j
## 1  11.04642 2770040 1500000 4270040      0.000000e+00 125000 2500000  11.73607
## 2  11.05880 2740012 1517357 4257370      2.328306e-10 126250 2470322  11.74602
## 3  11.07108 2710615 1534435 4245050      1.164153e-10 127500 2441274  11.75587
## 4  11.08325 2681831 1551238 4233069      0.000000e+00 128750 2412837  11.76563
## 5  11.09533 2653640 1567775 4221415      0.000000e+00 130000 2384992  11.77529
## 6  11.10730 2626025 1584052 4210077      1.164153e-10 131250 2357724  11.78486
##       v_j w_j     u_j expend_j-budget_j
## 1 2500012   0 2500012      0.000000e+00
## 2 2470333   0 2470333      1.164153e-10
## 3 2441285   0 2441285     -1.164153e-10
## 4 2412848   0 2412848     -1.164153e-10
## 5 2385004   0 2385004     -1.164153e-10
## 6 2357735   0 2357735      1.164153e-10
tail(mega_frame)
##      p1     p2      t1_i t2_i     psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j
## 1496  1 3.9875 0.3333333    1 0.04817813                    0    0    1
## 1497  1 3.9900 0.3333333    1 0.04847633                    0    0    1
## 1498  1 3.9925 0.3333333    1 0.04877490                    0    0    1
## 1499  1 3.9950 0.3333333    1 0.04907382                    0    0    1
## 1500  1 3.9975 0.3333333    1 0.04937307                    0    0    1
## 1501  1 4.0000 0.3333333    1 0.04967264                    0    0    1
##            z1       z2 doubleclearing     x1_i    x2_i        s1      s2
## 1496 -1834120 459967.3                -3824742 1709184 3127.4452 1499217
## 1497 -1818968 455881.7                -3811467 1705255 2501.5645 1499374
## 1498 -1803904 451823.1                -3798278 1701353 1875.8797 1499531
## 1499 -1788926 447791.3                -3785176 1697478 1250.3909 1499687
## 1500 -1774035 443786.1                -3772160 1693630  625.0977 1499844
## 1501 -1759230 439807.4                -3759230 1689807    0.0000 1500000
##              y1       y2 v_tilde_i v_i     w_i u_i expend_i-budget_i    x1_j
## 1496  -996872.6 499216.9       NaN NaN 3001561 NaN      4.656613e-10 1993750
## 1497  -997498.4 499373.8       NaN NaN 3001249 NaN     -1.396984e-09 1995000
## 1498  -998124.1 499530.6       NaN NaN 3000937 NaN      0.000000e+00 1996250
## 1499  -998749.6 499687.2       NaN NaN 3000625 NaN      9.313226e-10 1997500
## 1500  -999374.9 499843.7       NaN NaN 3000312 NaN      4.656613e-10 1998750
## 1501 -1000000.0 500000.0       NaN NaN 3000000 NaN     -4.656613e-10 2000000
##          x2_j v_tilde_j      v_j w_j      u_j expend_j-budget_j
## 1496 250000.6  14.50553 250015.1   0 250015.1      4.656613e-10
## 1497 250000.4  14.50615 250014.9   0 250014.9      0.000000e+00
## 1498 250000.2  14.50678 250014.7   0 250014.7     -4.656613e-10
## 1499 250000.1  14.50741 250014.6   0 250014.6      0.000000e+00
## 1500 250000.0  14.50803 250014.5   0 250014.5      4.656613e-10
## 1501 250000.0  14.50866 250014.5   0 250014.5      0.000000e+00
# mega_frame

# A useful piece of information can be, given the ts we chose, which followed from  how much they care in the linear part of utility about good 2 wrt good 1. i.e. the ratio t2/t1
mega_frame$t2_i[1]/mega_frame$t1_i[1]
## [1] 3

And we can see the whole economy when markets clear \(\pm 0.001N\) for the equilibria candidates below.

full_equilibria <- mega_frame[mega_frame$doubleclearing == "All clear!", ]  
print(full_equilibria, digits = 4)
##     p1 p2   t1_i t2_i  psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j        z1
## 301  1  1 0.3333    1 -0.3333                    0    0    1 8.996e-06
##             z2 doubleclearing  x1_i  x2_i    s1    s2 y1 y2 v_tilde_i    v_i
## 301 -8.996e-06     All clear! 5e+05 5e+05 1e+06 1e+06  0  0     13.12 666680
##       w_i     u_i expend_i-budget_i  x1_j  x2_j v_tilde_j   v_j w_j   u_j
## 301 3e+06 3666680                 0 5e+05 5e+05     13.12 5e+05   0 5e+05
##     expend_j-budget_j
## 301                 0
# print(full_equilibria, digits = 15)

Equilibrium (Best Guess) Description

In particular, our best guess of the equilirbium is as follows

loss=function(z1,z2){
  z1^2+z2^2}
loss_col=loss(full_equilibria$z1,full_equilibria$z2)
# print(loss_col)

# add loss column to mega frame
loss_mega=loss(mega_frame$z1,mega_frame$z2)
mega_frame=cbind(mega_frame,loss_mega)

# range of loss
range(mega_frame$loss_mega)
## [1] 1.618726e-10 1.499697e+32
full_equilibria = cbind(full_equilibria, loss_col)
# full_equilibria
full_equilibria[which.min(full_equilibria$loss_col),]
##     p1 p2      t1_i t2_i     psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j          z1
## 301  1  1 0.3333333    1 -0.3333333                    0    0    1 8.99646e-06
##               z2 doubleclearing  x1_i  x2_i    s1    s2 y1 y2 v_tilde_i
## 301 -8.99646e-06     All clear! 5e+05 5e+05 1e+06 1e+06  0  0  13.12236
##          v_i   w_i     u_i expend_i-budget_i  x1_j  x2_j v_tilde_j      v_j w_j
## 301 666679.8 3e+06 3666680                 0 5e+05 5e+05  13.12236 500013.1   0
##          u_j expend_j-budget_j     loss_col
## 301 500013.1                 0 1.618726e-10
# print(full_equilibria[which.min(full_equilibria$loss_col),][,c("p2","loss_col")], digits = 20)

p2_equilibrium=full_equilibria[which.min(full_equilibria$loss_col),]$p2
p2_equilibrium
## [1] 1

(If the market doesn’t clear (or if it clears more than once) our best equilibrium candidate is below.)

equilibrium_candidate = mega_frame


# add a column which directly targets loss when pointing out market clearing
low_loss_test=sapply(equilibrium_candidate$loss_mega, lowloss)
equilibrium_candidate =cbind(equilibrium_candidate, low_loss_test)

# equilibrium_candidate
equilibrium_candidate[which.min(equilibrium_candidate$loss_mega),]
##     p1 p2      t1_i t2_i     psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j          z1
## 301  1  1 0.3333333    1 -0.3333333                    0    0    1 8.99646e-06
##               z2 doubleclearing  x1_i  x2_i    s1    s2 y1 y2 v_tilde_i
## 301 -8.99646e-06     All clear! 5e+05 5e+05 1e+06 1e+06  0  0  13.12236
##          v_i   w_i     u_i expend_i-budget_i  x1_j  x2_j v_tilde_j      v_j w_j
## 301 666679.8 3e+06 3666680                 0 5e+05 5e+05  13.12236 500013.1   0
##          u_j expend_j-budget_j    loss_mega low_loss_test
## 301 500013.1                 0 1.618726e-10     Low loss!
equilibrium_candidate[equilibrium_candidate$low_loss_test == "Low loss!", ]
##     p1 p2      t1_i t2_i     psi1_i psi2_i=psi2_j=psi1_j t1_j t2_j          z1
## 301  1  1 0.3333333    1 -0.3333333                    0    0    1 8.99646e-06
##               z2 doubleclearing  x1_i  x2_i    s1    s2 y1 y2 v_tilde_i
## 301 -8.99646e-06     All clear! 5e+05 5e+05 1e+06 1e+06  0  0  13.12236
##          v_i   w_i     u_i expend_i-budget_i  x1_j  x2_j v_tilde_j      v_j w_j
## 301 666679.8 3e+06 3666680                 0 5e+05 5e+05  13.12236 500013.1   0
##          u_j expend_j-budget_j    loss_mega low_loss_test
## 301 500013.1                 0 1.618726e-10     Low loss!

Behaviour at around lowest loss

# can nicely print some relevant facts to observe low loss behaviour by uncommenting line below
# print(equilibrium_candidate[,c(1,2,9,10,11,30,31)], digits = 5)
# to get our equilibria candidates we split the prices into sections (5 below but variable) and add to a vector, the minimum loss of that candidate (this is using local uniqueness but can run into problems with certain functions. an extension of the work here is to generalise the way this is done)
# chunk <- function(x,n) split(x, cut(seq_along(x), n, labels = FALSE)) #splitting any x into n groups
num_chunks=5
size_of_chunks= ceiling((num_of_steps-1)/num_chunks)
# size_of_chunks

# multiple_equilibria
# set equilibirum to be this instead
min_local_loss = function(df){
 df[which.min(df$loss_mega),]$p2
}
min_local_loss(equilibrium_candidate)
## [1] 1
# dummy
group_dummy=c()
for (i in 1:num_chunks) {
  group_dummy=c(group_dummy,rep(i,size_of_chunks))
}
# add whatever extra to the last group (usually one)
group_dummy=c(group_dummy,rep(num_chunks,num_of_steps-length(group_dummy)))
equilibrium_candidate=cbind(equilibrium_candidate,group_dummy)

# num_of_steps-length(group_dummy) # should be zero
# print(equilibrium_candidate[,c(2,11,30:32)], digits = 5)




# candidate lacking clearing... 
p2_equilibrium_candidate=min_local_loss(equilibrium_candidate)
p2_equilibrium_candidate
## [1] 1
# key multiple equilibria candidates
price_with_min_loss_per_group=daply (mega_frame, .(group_dummy), min_local_loss)# daply: Split data frame, apply function, and return results in an array.
# the p2_equilirbium/a below is the intersection of those eqa with low loss and that minimise losss locally (in their groups)
p2_equilibrium=intersect(price_with_min_loss_per_group, equilibrium_candidate[equilibrium_candidate$low_loss_test == "Low loss!", ]$p2)

We thus have the following equilibria \(\bar{p}_2=\)

p2_equilibrium
## [1] 1

Appendix

Delta bar

We supposed that \(\bar{\delta} = e_2/\bar{p}_2\).

Check that it’s respected.

\[ \chi_1^i(p) = \frac{p_2}{p_1(t_2^i+\psi_2^i) -p_2(t_1^i + \psi_1^i)} \]

Which gives zero derivative with respect to profit.

\[ \chi_2^i(p) = \frac{p_1}{p_2}\left(\omega_1^i - \frac{p_2}{p_1(t_2^i+\psi_2^i) - p_2(t_1^i + \psi_1^i)}\right) + \omega_2^i + \frac{\pi(p)}{Np_2} \]

Which indeed gives \(\frac{1}{Np_2}\), so when aggregated it gives \(\bar{\delta} = e_2/\bar{p}_2\).

At equilibrium:

bar_delta=c(0,1/(N*full_equilibria[which.min(full_equilibria$loss_col),]$p2)) #this is delta for each individual!
bar_delta
## [1] 0e+00 1e-06

Ensuring \(F\) works as expected at equilibria

ginv_tests =function(){
   message("Summary of Key parameters")
print(as.data.frame(matrix(c("t1_i","t2_i","psi1_i","psi2_i",t1,t2,"variable",psi2_i,"t1_j","t2_j","psi1_j","psi2_j",t1_j,t2_j,psi1_j,psi2_j),nrow=4,byrow = TRUE), row.names = NULL, optional = FALSE, col.names=c("1,22")))


# psi1=nleqslv(psi1start, fnpsi_minuspsi, control=list(btol=.01))$x     if we need to check psi1
  
for (i in 1:length(p2_equilibrium)){
  p2=p2_equilibrium[i]
psi1=equilibrium_candidate[equilibrium_candidate$p2 == p2_equilibrium[i], ]$psi1_i

  print(" c(p1,p2,t1,t2,psi1,psi2_i,p1,p2,t1_j,t2_j,psi1_j,psi2_j)")
  print(c(p1,p2,t1,t2,psi1,psi2_i,p1,p2,t1_j,t2_j,psi1_j,psi2_j))
  Jx=J_x(p1,p2,t1,t2,psi1,psi2_i)+J_x(p1,p2,t1_j,t2_j,psi1_j,psi2_j)#J_x for j which is what i best responds to when calculating -Jz leading to F
  Js=J_s(1,p2)
  minus_Jz=Js-Jx
  F_inv=F(Js,Jx)
  
  
  message("----------------------------------------")
  message("") #this adds an empty line
  message(paste0(paste0("For eqilibrium ",i),"------------------------"))
  # print("----------------------------------------", quote = FALSE) #quote = FALSE)gets rid of the " " 
  cat(paste0("We have p2 = ",p2))
  
      message("") 
  cat("J_s=")
  print(Js)
  
  cat("Jx=")
  print(Jx)
  cat("-Jz=")
  print(minus_Jz)
  
      message("")
  cat("F_inv=")
  print(F_inv)
      message("")

  # AA^gA=A for A^g the generalised inverse of A
  # so the code below should give -Jz
  message("-Jz %*% F_inv %*% -Jz = (should be -Jz)")
  print(minus_Jz %*% F_inv %*% minus_Jz)
  message("-Jz %*% F_inv %*% -Jz should be -Jz giving zeros below")
  print(minus_Jz %*% F_inv %*% minus_Jz - minus_Jz)
  print(closetozero(norm(minus_Jz %*% F_inv %*% minus_Jz - minus_Jz ),0))

  # You can test your own F matrices guesses to check they are working by changing Ftest below
  Ftest=matrix(c(1/3,0,0,0),nrow=2)
  
  message("You picked Ftest=")
  print(Ftest)
  message("-Jz %*% Ftest %*% -Jz = (should be -Jz)")
  print(minus_Jz %*% Ftest %*% minus_Jz)
  message("-Jz %*% Ftest %*% -Jz should be -Jz giving zeros below")
  closetozero(norm(minus_Jz %*% Ftest %*% minus_Jz - minus_Jz),0)
  
  message("-------psi tests---------")

  
  message("check def of psis: psis=(Js F)^t grad_externalities")
  print(t(t( Js  %*% F_inv) %*% t(grad_externalities2(p1,p2))))
  psi1_at_eqm=equilibrium_candidate[equilibrium_candidate$p2 == p2, ]$psi1
  psis_test=as.data.frame(matrix(c(psi1_at_eqm,psi2_i),nrow=1,byrow = TRUE))
  colnames(psis_test)=c("psi1_i","psi2_i")
  print(psis_test)
  
  message("If test above is true the line below should tend to 1")
  print({t( Js  %*% F_inv) %*% t(grad_externalities2(p1,p2))}[1]/psi1_at_eqm)
  print({t( Js  %*% F_inv) %*% t(grad_externalities2(p1,p2))}[1]/psi1_at_eqm <1.01 &{t( Js  %*% F_inv) %*% t(grad_externalities2(p1,p2))}[1]/psi1_at_eqm >0.99) #Test true or false
  
  message("-------delta tests---------")
  
  bar_delta=c(0,1/p2)
  cat("bar_delta=")
  print(bar_delta)
  
  message("check Fdelta_bar=zeros")
  print(F_inv %*% bar_delta)
}
}
ginv_tests()
Summary of Key parameters
                 V1   V2       V3     V4
1              t1_i t2_i   psi1_i psi2_i
2 0.333333333333333    1 variable      0
3              t1_j t2_j   psi1_j psi2_j
4                 0    1        0      0
[1] " c(p1,p2,t1,t2,psi1,psi2_i,p1,p2,t1_j,t2_j,psi1_j,psi2_j)"
 [1]  1.0000000  1.0000000  0.3333333  1.0000000 -0.3333333  0.0000000
 [7]  1.0000000  1.0000000  0.0000000  1.0000000  0.0000000  0.0000000
----------------------------------------
For eqilibrium 1------------------------
We have p2 = 1
J_s=       [,1]   [,2]
[1,]  5e+05 -5e+05
[2,] -5e+05  5e+05
Jx=       [,1]   [,2]
[1,] -1e+06  1e+06
[2,]  1e+06 -1e+06
-Jz=         [,1]     [,2]
[1,]  1500000 -1500000
[2,] -1500000  1500000
F_inv=              [,1]          [,2]
[1,]  6.666667e-07  0.000000e+00
[2,] -7.486784e-23 -2.245194e-34
-Jz %*% F_inv %*% -Jz = (should be -Jz)
         [,1]     [,2]
[1,]  1500000 -1500000
[2,] -1500000  1500000
-Jz %*% F_inv %*% -Jz should be -Jz giving zeros below
              [,1]          [,2]
[1,]  9.313226e-10 -9.313226e-10
[2,] -9.313226e-10  9.313226e-10
[1] "All clear!"
You picked Ftest=
          [,1] [,2]
[1,] 0.3333333    0
[2,] 0.0000000    0
-Jz %*% Ftest %*% -Jz = (should be -Jz)
         [,1]     [,2]
[1,]  7.5e+11 -7.5e+11
[2,] -7.5e+11  7.5e+11
-Jz %*% Ftest %*% -Jz should be -Jz giving zeros below
-------psi tests---------
check def of psis: psis=(Js F)^t grad_externalities
           [,1]          [,2]
[1,] -0.3333333 -1.122597e-28
      psi1_i psi2_i
1 -0.3333333      0
If test above is true the line below should tend to 1
[1] 1
[1] TRUE
-------delta tests---------
bar_delta=[1] 0 1
check Fdelta_bar=zeros
              [,1]
[1,]  0.000000e+00
[2,] -2.245194e-34

Graph of \(x_1\) and of \(x_2\) as a function of the prices with \(p_1=1\)

ggplot(data = mega_frame, mapping = aes(x = p2, y = x1_i, color="x1_i"))+
  geom_point(size = 0.5, alpha=0.05)+geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "purple") + annotate("text", x=p2_equilibrium-0.05, y=legend_position(min(mega_frame$x2_i),max(mega_frame$x2_i)), label="at equilibrium", angle=90, color = "purple",size = 3.5) +
  geom_point(data = mega_frame, mapping = aes(x = p2, y = x2_i, color="x2_i"), alpha=0.05)+
  geom_point(data = mega_frame, mapping = aes(x = p2, y = x1_j, color="x1_j"), alpha=0.05)+
  geom_point(data = mega_frame, mapping = aes(x = p2, y = x2_j, color="x2_j"), alpha=0.05)+
  labs(x = TeX("$p_2$"), y = TeX("$x^k_\U2113$")) + theme_light() +ylim(0,N)
## Warning: Removed 1055 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1148 rows containing missing values (geom_point).
## Warning: Removed 800 rows containing missing values (geom_point).
## Warning: Removed 115 rows containing missing values (geom_point).

# +  annotate_eqm_grey('x2_i') + grey_line

Graph of \(\psi^i_1\) as a function of the prices with \(p_1=1\)

ggplot(data = mega_frame, mapping = aes(x = p2, y = psi1_i)) +
  geom_point(size = 0.5)+geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "purple") +annotate("text", x=p2_equilibrium-0.05, y=legend_position(min(mega_frame$psi1),max(mega_frame$psi1)), label="at equilibrium", angle=90, color = "purple",size = 3.5)+ labs(x = TeX("$p_2$"), y = TeX("$\\psi^i_1$"))

# + annotate_eqm_grey('psi1') + grey_line

Graph of \(s_1\) and \(s_2\) as a function of the prices with \(p_1=1\)

ggplot(data = mega_frame, mapping = aes(x = p2, y = s1, color="s1")) +
  geom_point(size = 0.5)+geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "purple") +annotate("text", x=p2_equilibrium-0.05, y=legend_position(min(mega_frame$s1),max(mega_frame$s1)), label="at equilibrium", angle=90, color = "purple",size = 3.5) + geom_point(data = mega_frame, mapping = aes(x = p2, y = s2, color="s2"))+ labs(x = TeX("$p_2$"), y = TeX("$s_\U2113$$")) + theme_light()

# + annotate_eqm_grey('s1') + grey_line

Graph of \(z_1\) and of \(z_2\) as a function of the prices with \(p_1=1\)

ggplot(data = mega_frame, mapping = aes(x = p2, y = z1, color="z1")) +
  # geom_line(y=z2, color="red") +
  geom_point(size = 0.5)+geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "purple") +annotate("text", x=p2_equilibrium-0.05, y=legend_position(min(mega_frame$z2),max(mega_frame$z2)), label="at equilibrium", angle=90, color = "purple",size = 3.5)+ geom_point(data = mega_frame, mapping = aes(x = p2, y = z2, color="z2"))+ labs(x = TeX("$p_2$"), y = TeX("$z_\U2113$")) + theme_light()+ylim(-N,N)
## Warning: Removed 1107 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1016 rows containing missing values (geom_point).

# + annotate_eqm_grey('z2') + grey_line

[Zoom]

p2_equilibrium2=p2_equilibrium[1]
ggplot(data = mega_frame, mapping = aes(x = p2, y = z1, color="z1")) + geom_point(size = 0.5) +annotate("text", x=p2_equilibrium2-0.05, y=legend_position(min(mega_frame$z2),max(mega_frame$z2)), label="at equilibrium", angle=90, color = "purple",size = 3.5)+ geom_point(data = mega_frame, mapping = aes(x = p2, y = z2, color="z2"))+ labs(x = TeX("$p_2$"), y = TeX("$z_\U2113$")) + theme_light() + xlim(p2_equilibrium2-0.1, p2_equilibrium2+0.1)+ylim(-N/10,N/10) +geom_vline(xintercept = p2_equilibrium2, linetype= "dashed", color = "purple")

[Zoom2]

p2_equilibrium1=p2_equilibrium[2]
ggplot(data = mega_frame, mapping = aes(x = p2, y = z1, color="z1")) + geom_point(size = 0.5) +annotate("text", x=p2_equilibrium1-0.05, y=legend_position(min(mega_frame$z2),max(mega_frame$z2)), label="at equilibrium", angle=90, color = "purple",size = 3.5)+ geom_point(data = mega_frame, mapping = aes(x = p2, y = z2, color="z2"))+ labs(x = TeX("$p_2$"), y = TeX("$z_\U2113$")) + theme_light() + xlim(p2_equilibrium1-0.1, p2_equilibrium1+0.1)+ylim(-N/100,N/100) +geom_vline(xintercept = p2_equilibrium1, linetype= "dashed", color = "purple")

From this point on, \(i\) will have an light green utility curve and \(j\) a light blue one.

Graph of non linear consumption utility \(\tilde{v}^k_1(x_1)=\ln(x_1)\) as a function of the prices with \(p_1=1\).

ggplot(data = mega_frame, mapping = aes(x = p2, y = v_tilde_i, color="v_tilde_i")) +
  geom_point(size = 0.5,color = "lightgreen")+
  geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue")+ theme_bw() +
  labs(x = TeX("$p_2$"), y = TeX("$\\tilde{v}^k$")) +
  annotate("text", x=p2_equilibrium-0.05, y=legend_position(min(mega_frame$v_tilde_i),max(mega_frame$v_tilde_i)), label="at equilibrium", angle=90, color = "darkblue",size = 3.5)+ geom_point(data = mega_frame, mapping = aes(x = p2, y = v_tilde_j, color="v_tilde_j"), color = "lightblue")
## Warning: Removed 401 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).

# + annotate_eqm_grey('v_tilde') + grey_line
# i likes good 1 less because of the supply externality and hence the curve v_tilde of utility derived from x1 is below that for j

Graph of consumption utility \(v^k(x) = \tilde{v}^k_1(x_1) + t^k \cdot x\) as a function of the prices with \(p_1=1\).

ggplot(data = mega_frame, mapping = aes(x = p2, y = v_i, color="v_i")) +
  geom_point(size = 0.5,color = "lightgreen")+geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue")+ theme_bw() +annotate("text", x=p2_equilibrium-0.05,  y=legend_position(min(mega_frame$v_i),max(mega_frame$v_i)), label="at equilibrium", angle=90, color = "darkblue",size = 3.5) +
      labs(x = TeX("$p_2$"), y = TeX("$v^k$")) +
  geom_point(data = mega_frame, mapping = aes(x = p2, y = v_j, color="v_j"), color = "lightblue")
## Warning: Removed 401 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).

# + annotate_eqm_grey('v') + grey_line

Graph of externality-concerns utility \(w^k(s_1,s_2)\) as a function of the prices with \(p_1=1\).

ggplot(data = mega_frame, mapping = aes(x = p2, y = w_i,color="w_i")) +
  geom_point(size = 0.5,color = "lightgreen")+geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue")+ theme_bw() +annotate("text", x=p2_equilibrium-0.05,  y=legend_position(min(min(mega_frame$w_i),min(mega_frame$w_j)),max(max(mega_frame$w_i),max(mega_frame$w_j))), label="at equilibrium", angle=90, color = "darkblue",size = 3.5) +
    geom_point(data = mega_frame, mapping = aes(x = p2, y = w_j, color="w_j"), color = "lightblue")+
    labs(x = TeX("$p_2$"), y = TeX("$w^k$"))

# + annotate_eqm_grey('w') + grey_line

Graph of full utility \(u^k(x,s)=v^k(x)+w^k(s)\) as a function of the prices with \(p_1=1\).

color_vari=c(1,1)
ggplot(data = mega_frame, mapping = aes(x = p2, y = u_i, color=color_vari)) +
  geom_point(size = 0.5,color =  "lightgreen")+
  geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue")+ theme_bw() +
  annotate("text", x=p2_equilibrium-0.05,
           y=legend_position(min(mega_frame$u_i),max(mega_frame$u_i)),
           label="at equilibrium", angle=90,
           color = "darkblue",size = 3.5) +
  geom_point(data = mega_frame, mapping = aes(x = p2, y = u_j), color = "lightblue")+
    labs(x = TeX("$p_2$"), y = TeX("$u^k$"), color="Colour Key")
## Warning: Removed 401 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_text).

  # scale_color_manual(labels = c("T999", "T888"), values = c("blue", "red"))
  # annotate_eqm_grey('u') + grey_line

Graph of loss as a function of the prices with \(p_1=1\).

color_vari=c(1,1)
ggplot(data = mega_frame, mapping = aes(x = p2, y = loss_mega, color=color_vari)) +
  geom_point(size = 0.5,color =  "purple")+
  geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue",alpha=0.2)+ theme_bw() +
  annotate("text", x=p2_equilibrium-0.05,
           y=legend_position(min(mega_frame$loss_mega),N*100000),
           label="at equilibrium", angle=90,
           color = "darkblue",size = 3.5,alpha=0.2) +
    labs(x = TeX("$p_2$"), y = "loss", color="Colour Key") + ylim(min(mega_frame$loss_mega),N*100000)
## Warning: Removed 1412 rows containing missing values (geom_point).

  # scale_color_manual(labels = c("T999", "T888"), values = c("blue", "red"))
  # annotate_eqm_grey('u') + grey_line

[Zoom] Graph of loss as a function of the prices with \(p_1=1\).

color_vari=c(1,1)
ggplot(data = mega_frame, mapping = aes(x = p2, y = loss_mega, color=color_vari)) +
  geom_point(size = 2,color =  "purple")+
  geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue",alpha=0.2)+ theme_bw() +
  annotate("text", x=p2_equilibrium-0.05,
           y=legend_position(min(mega_frame$loss_mega),N*1000),
           label="at equilibrium", angle=90,
           color = "darkblue",size = 3.5,alpha=0.2) +
    labs(x = TeX("$p_2$"), y = "loss", color="Colour Key") + ylim(min(mega_frame$loss_mega),N*1000)
## Warning: Removed 1492 rows containing missing values (geom_point).

  # scale_color_manual(labels = c("T999", "T888"), values = c("blue", "red"))
  # annotate_eqm_grey('u') + grey_line

[ZoomOut] Graph of loss as a function of the prices with \(p_1=1\).

color_vari=c(1,1)
ggplot(data = mega_frame, mapping = aes(x = p2, y = loss_mega, color=color_vari)) +
  geom_point(size = 2,color =  "purple")+
  geom_vline(xintercept = p2_equilibrium, linetype= "dashed", color = "darkblue",alpha=0.2)+ theme_bw() +
  annotate("text", x=p2_equilibrium-0.05,
           y=legend_position(min(mega_frame$loss_mega),N*1000),
           label="at equilibrium", angle=90,
           color = "darkblue",size = 3.5,alpha=0.2) +
    labs(x = TeX("$p_2$"), y = "loss", color="Colour Key") + ylim(min(mega_frame$loss_mega),max(mega_frame$loss_mega)/(N*1e+15))+xlim (p2_equilibrium[1]-0.1,p2_equilibrium[1]+0.1)
## Warning: Removed 1420 rows containing missing values (geom_point).

  # scale_color_manual(labels = c("T999", "T888"), values = c("blue", "red"))
  # annotate_eqm_grey('u') + grey_line