Solve in Maxima:

Of course, we should segregate wastes, but the problem is that the industry produces too many plastic packing. Why? Because the demand of the consumer-related industries is (and will be) growing fast. The production of plastic packing should be constrained.

Let's analyse the problem using the simple partial equilibrium model.

We have three collective actors:

* the packing industry - the supplier of plastic packing

* consumer-related industries - demand side of the plastic packing market.

* consumers (society) who are the recipients of externalities of the production of the packing industry. Externalities can be positive or negative. Packing protects consumer goods - positive. Littering the earth - negative. After achieving a saturation point, the production of plastic packing generates negative externalities.

Let's assume that the partial equilibrium model is an approximation of market conditions:

demand: q = 200*N - 3pd = supply: q = -100 + 2ps equilibrium: pd = ps

Production of plastic packing generates the externalities EX which can be approximated by a function: EX = 120*q - 0.04*q^2

b) Society stops to be passive and start to boycott the plastic packing. The boycott causes a decrease of demand by the factor b.

demand: q = b*(200*N - 3p_d ) = .. [b > 0 and b <1]

supply: q = -100 + 2p_s

equilibrium: p_d = p_s

How the boycott influence social welfare when society wants to reduce the externality effect to zero EX = 0 (for q>0)? Compare with the benchmark situation. Analyze the obtained results. Determine the impact of the boycott on this market. (11p)

[Hint: present: equilibrium condition (q*, pd* = ps*) , consumer surplus (CS), producer surplus (PS), externalities EX, Social Welfare SW = CS + SP + EX as functions of b and graph the functions. Modify the Maxima codes].

kill(all)$ load(draw)$ assume(p_s>0,p_d>0,q>0, q_star>0)$ /* Demand function*/ eq1:q = b*(1000 - 5*p_d)$ solve([eq1],[p_d])$ p_d:rhs(%[1])$ /* Inverse demand function */ /* Suplly function */ eq2: q = (-250 + 2*p_s)$ solve([eq2],[p_s])$ p_s:rhs(%[1])$ /* Inverse suplly function */ /* Equlibrium */ eq3: p_d = p_s$ solve([eq3],[q])$ q_star:rhs(%[1])$ q:ev(''q_star)$ p_s_star:ev(''p_s)$ p_d_star:ev(''p_d)$ kill(q)$ CS: ev(integrate(''p_d, q, 0, ''q_star) - integrate(''p_d_star, q, 0, ''q_star))$ PS: ev(integrate(''p_s_star, q, 0, ''q_star) - integrate (''p_s, q, 0, ''q_star))$ EX : - 0.5*''q_star^2 + 30*q_star $ SW: ev(''CS + ''PS + EX)$ draw2d( grid= true, xlabel = "b", ylabel = "PS, CS, EX, SW", color = grey, color = red, key = " CS ", explicit(''CS,b,0,1), color = light-green, key = " PS ", explicit(''PS,b,0,1), color = light-blue, key = " SW ", explicit(''SW,b,0,1), color = brown, key = " EX ", explicit(''EX,b,0,1) ) $

kill(all)$ load(draw)$ assume(p_s>0,p_d>0,q>0,q>0,q_star>0)$ /* Demand function*/ eq1:q = b*(1000 - 5*p_d)$ solve([eq1],[p_d])$ p_d:rhs(%[1])$ /* Inverse demand function */ /* Suplly function */ eq2:q = (-250 + 2*p_s)$ solve([eq2],[p_s])$ p_s:rhs(%[1])$ /* Inverse suplly function */ /* Equlibrium */ eq3: p_d = p_s$ solve([eq3],[q])$ q_star:rhs(%[1])$ q:ev(''q_star)$ p_s_star:ev(''p_s)$ p_d_star:ev(''p_d)$ kill(q)$ EX : - 0.5*''q_star^2 + 30*q_star $ eq1: ''EX = 0; solve(eq1,b); b_0: float(rhs(%[1])); print("*****results***************" )$ print("EX = 0 for b = ", float(b_0))$ kill(all)$ load(draw)$ assume(p_s>0,p_d>0,q>0,q>0,q_star>0)$ b: 0.267; /* insert b */ /* Demand function*/ eq1:q = b*(1000 - 5*p_d)$ solve([eq1],[p_d])$ p_d:rhs(%[1])$ /* Inverse demand function */ /* Suplly function */ eq2:q = (-250 + 2*p_s)$ solve([eq2],[p_s])$ p_s:rhs(%[1])$ /* Inverse suplly function */ /* Equlibrium */ eq3: p_d = p_s$ solve([eq3],[q])$ q_star:rhs(%[1])$ q:ev(''q_star)$ p_s_star:ev(''p_s)$ p_d_star:ev(''p_d)$ kill(q)$ CS: ev(integrate(''p_d, q, 0, ''q_star) - integrate(''p_d_star, q, 0, ''q_star)); PS: ev(integrate(''p_s_star, q, 0, ''q_star) - integrate (''p_s, q, 0, ''q_star)); EX: 500*''q_star - 0.01*''q_star^2; SW: ev(''CS + ''PS + EX); print("*****results***************" )$ print("q = ", float(q_star) )$ print("p_d = ", float(p_d_star) )$ print("p_s = ", float(p_s_star) )$ print("PS = ", float(PS) )$ print("CS = ", float(CS) )$ print("EX = ", float(EX) )$ print("SW = ", float(SW))$