# ECON 331 97/1
# ASSIGNMENT #9
# 
# QUESTION 1
# 
# The general form of the slutsky equation is:  Diff(x,Px) =
# -x(Diff(x,B)) + Diff(x^c,Px) , where x^c is the compensated demand
# function.
# 
# What follows is certainly not the most efficient way to verify the
# Slutsky equations for a specific utility function.  The purpose here
# is to follow the procedure outlined in Chiang, pp. 404-408.
# 
> restart;
# 
# Set up the Lagrangian:
# 
> U:=x*y:
> L:=U+lambda*(B-Px*x-Py*y);

                 L := x y + lambda (B - Px x - Py y)

# 
# First-order conditions:
# 
> focs:=(diff(L,lambda),diff(L,x),diff(L,y));

        focs := B - Px x - Py y, y - lambda Px, x - lambda Py

# 
# Comparative Statics:
# 
> with(linalg):
> H:=hessian(L,[lambda,x,y]);
Warning, new definition for norm
Warning, new definition for trace

                            [ 0     -Px    -Py]
                            [                 ]
                       H := [-Px     0      1 ]
                            [                 ]
                            [-Py     1      0 ]

# 
# Solution vector for differentiating with respect to B:
# 
> a:=map(x->-x,jacobian([focs[1],focs[2],focs[3]],[B]));

                                   [-1]
                                   [  ]
                              a := [ 0]
                                   [  ]
                                   [ 0]

# 
# Comparative statics results for B:
# 
> matrix(3,1,[Diff(lambda,B),Diff(x,B),Diff(y,B)])=linsolve(H,a);

                      [d        ]   [      1  ]
                      [-- lambda]   [1/2 -----]
                      [dB       ]   [    Py Px]
                      [         ]   [         ]
                      [  d      ]   [     1   ]
                      [  -- x   ] = [1/2 ---- ]
                      [  dB     ]   [     Px  ]
                      [         ]   [         ]
                      [  d      ]   [     1   ]
                      [  -- y   ]   [1/2 ---- ]
                      [  dB     ]   [     Py  ]

> sols_B:=linsolve(H,a):
# 
# Solution vector for differentiating with respect to P[x]:
# 
> c:=map(x->-x,jacobian([focs[1],focs[2],focs[3]],[Px]));

                                 [  x   ]
                                 [      ]
                            c := [lambda]
                                 [      ]
                                 [  0   ]

# 
# Comparative static results for P[x]:
# 
> matrix(3,1,[Diff(lambda,Px),Diff(x,Px),Diff(y,Px)])=linsolve(H,c);

                 [ d        ]   [      lambda Py + x]
                 [--- lambda]   [- 1/2 -------------]
                 [dPx       ]   [          Py Px    ]
                 [          ]   [                   ]
                 [   d      ]   [      lambda Py + x]
                 [  --- x   ] = [- 1/2 -------------]
                 [  dPx     ]   [           Px      ]
                 [          ]   [                   ]
                 [   d      ]   [     lambda Py - x ]
                 [  --- y   ]   [ 1/2 ------------- ]
                 [  dPx     ]   [          Py       ]

> sols_Px:=linsolve(H,c):
# Solution vector for differentiating with respect to P[x] with
# diff(B,P[x])=x: 
> d:=map(x->-x,matrix(3,1,[0,-lambda,0]));

                                 [  0   ]
                                 [      ]
                            d := [lambda]
                                 [      ]
                                 [  0   ]

# 
# Comparative static results for the `compensated` demand functions:
# 
> matrix(3,1,[Diff(lambda^c,Px),Diff(x^c,Px),Diff(y^c,Px)])=linsolve(H,d
> );

                  [ d        c]
                  [--- lambda ]   [       lambda  ]
                  [dPx        ]   [ - 1/2 ------  ]
                  [           ]   [         Px    ]
                  [   d   c   ]   [               ]
                  [  --- x    ] = [      Py lambda]
                  [  dPx      ]   [- 1/2 ---------]
                  [           ]   [         Px    ]
                  [   d   c   ]   [               ]
                  [  --- y    ]   [  1/2 lambda   ]
                  [  dPx      ]

> csols_Px:=linsolve(H,d):
# 
# The components of the two Slutsky equations are:
# 
> Inc_x:=-x*sols_B[2,1];

                                         x
                         Inc_x := - 1/2 ----
                                         Px

> Sub_x:=csols_Px[2,1];

                                      Py lambda
                       Sub_x := - 1/2 ---------
                                         Px

> Inc_y:=-x*sols_B[3,1];

                                         x
                         Inc_y := - 1/2 ----
                                         Py

> Sub_y:=csols_Px[3,1];

                         Sub_y := 1/2 lambda

# 
# We can sign these effects by substituting the actual functions which
# are obtained from the first-order conditions:
# 
> sols:=solve({diff(L,lambda),diff(L,x),diff(L,y)},{lambda,x,y});

                               B             B             B
       sols := {lambda = 1/2 -----, y = 1/2 ----, x = 1/2 ----}
                             Px Py           Py            Px

# 
# Substituting these solutions into the respective substitution and
# income effects, it is clear that the income and substitution effects
# for x, and the income effect for y are negative, while the
# substitution effect for y is positive:
# 
> 'Inc_x'=Inc_x;
> 'Sub_x'=Sub_x;
> 'Inc_y'=Inc_y;
> 'Sub_y'=Sub_y;
# 

                                         x
                          Inc_x = - 1/2 ----
                                         Px


                                     Py lambda
                       Sub_x = - 1/2 ---------
                                        Px


                                         x
                          Inc_y = - 1/2 ----
                                         Py


                          Sub_y = 1/2 lambda

# 
# 
# QUESTION 2
# 
# a)
> restart;
# 
# The limit-output model:
# 
> p:=120-(q1+q2):
> c1:=75+35*q1:
> c2:=100+40*q2:
> pi1:=p*q1-c1;

                pi1 := (120 - q1 - q2) q1 - 75 - 35 q1

> pi2:=p*q2-c2;

               pi2 := (120 - q1 - q2) q2 - 100 - 40 q2

# 
# The best response of firm 2 to the ouput produced by firm 1:
# 
> q2:=solve(diff(pi2,q2),q2);

                          q2 := 40 - 1/2 q1

# 
# Firm 1 chooses output (lim_q1) so that pi[2]=0 :
> lim_q1:=solve(pi2,q1);

                          lim_q1 := 60, 100

# 
# Choose the value for lim_q1 which maximises pi[1]:
> q1:=lim_q1[1]:
> q2:=0:
> 'pi1(lim_q1=60)'=pi1;

                       pi1(lim_q1 = 60) = 1425

> 'p(lim_q1=60)'=p;

                         p(lim_q1 = 60) = 60

> q1:=lim_q1[2]:
> 'pi1(lim_q1=100)'=pi1;
# 

                      pi1(lim_q1 = 100) = -1575

# So, lim_q1 = 60.
# 
# Set up functions for plot:
# 
> p:=120-Q:
> q1:=lim_q1[1]:
# 
# Firm 2 faces the 'residual' demand (D2 in plot):
# 
> q2:=Q-lim_q1[1];

                             q2 := Q - 60

> ac2:=c2/q2;

                                -2300 + 40 Q
                         ac2 := ------------
                                   Q - 60

> mr2:=60-2*q2;

                           mr2 := 180 - 2 Q

> mc2:=40:
# 
> with(plots):
> Z:=plot({p,mr2,ac2,mc2},Q=60.1..100, P=0..70,title=`Limit Output
> Equilibrium`):
> U:=textplot([100,46,'AC2']):
> V:=textplot([100,40,'MC2']):
> W:=textplot([100,18,'D2']):
> J:=textplot([90,5,'MR2']):
> display([Z,U,V,W,J]);
# 
# 

# 
# b)
> restart;
# 
# The Cournot Duopoly:
# 
> p:=120-(q1+q2):
> c2:=100+40*q2:
> c1:=75+35*q1:
> pi1:=p*q1-c1:
> pi2:=p*q2-c2:
# 
# The Cournot quantities:
# 
> cournot_q:=solve({diff(pi1,q1),diff(pi2,q2)},{q1,q2});

                   cournot_q := {q2 = 25, q1 = 30}

# Cournot price:
# 
> cournot_p:=subs(cournot_q[1],cournot_q[2],p);

                           cournot_p := 65

# Cournot profits:
# 
> profit[1]:=subs(cournot_q[1],cournot_q[2],pi1);

                           profit[1] := 825

> profit[2]:=subs(cournot_q[1],cournot_q[2],pi2);

                           profit[2] := 525

> with(plots):
> Z:=implicitplot({diff(pi1,q1),diff(pi2,q2),pi1=825,pi2=525},q1=0..100,
> q2=0..100,title=`Cournot Equilibrium`):
> U:=textplot([5,80,'R1']):
> V:=textplot([32,40,'Isoprofit_2']):
> W:=textplot([5,40,'R2']):
> J:=textplot([45,25,'Isoprofit_1']):
> display([Z,U,V,W,J]);
# 

# 
# c)
# 
# The Stackelberg Model:
# 
# Firm 2 best response function:
# 
> q2:=solve(diff(pi2,q2),q2);

                          q2 := 40 - 1/2 q1

# Firm 1 maximises profit subject to the best response from firm 2:
# 
# Substitute the best response function into pi[1]:
> pi1;

                    (80 - 1/2 q1) q1 - 75 - 35 q1

# 
# Solve for the Stackelberg q1 from the first-order condition:
# 
> q1:=solve(diff(pi1,q1),q1);

                               q1 := 45

# Stackelberg q2:
# 
> q2;

                                 35/2

# Stackelberg profits:
# 
> pi1;

                                1875/2

> pi2;

                                825/4

# Stackelberg price:
# 
> p;

                                115/2

> restart;
> p:=120-(q1+q2):
> c1:=75+35*q1:
> c2:=100+40*q2:
> pi1:=p*q1-c1:
> pi2:=p*q2-c2:
> with(plots):
> Z:=implicitplot({pi1=1875/2,diff(pi2,q2)=0},q1=0..80,q2=0..40,title=`S
> tackelberg Equilibrium`):
> U:=textplot([25,30,'R2']):
> V:=textplot([20,10,'Isoprofit_1']):
> display([Z,U,V]);
# 

