Generalizing to N variables

Question

I have a function of the following type where $p(i)$ is a variable on which the agent can optimize and $p(-i)$ a vector of variables that he has to take as parameters.

This is a demand function and from this function, one can compute a profit function.

My goal is to find the optimal solutions $p_1^*,...,p_N^*$ by computing the first-order conditions of each profit function with respect to $p(i)$ and setting each of them to 0.

Following the first answer of @user42582, I modified my code as follows:

(* Demand function *)
(* Takes a vector p1,...,pn and optional scalar arguments v and g and
returns a vector q1, ..., qn *)
(* Rn to Rn *)
q[p_, v_: V, g_: \[Gamma]] := Module[
  {meanP = Mean[p]},
  v - p - g (p - meanP)]
(* Profit function *)
(* Takes a vector p1,...,pn and optional scalar arguments cost, v and \
g and returns a vector profit1,...,profitn *)
(* Rn to Rn *)
profit[p_, cost_: c, v_: V, g_: \[Gamma]] := (p - cost) * q[p, v, g]

(* Getting the Jacobian matrix *)
(* Differentiating a Rn to Rn function wrt to X => N x N matrix *)
dvector[p_, cost_: c, v_: V, g_: \[Gamma]] := 
 D[profit[p, cost, v, g], {p}]
(* Displaying the Jacobian on an example *)
MatrixForm[dvector[{x, y, z}]]

(* Displaying only the element 2,2 of the Jacobian matrix *)
dvector[{x, y, z}][[2, 2]]

(* Solving for the prices *)
solution[p_, cost_: c, v_: V, g_: \[Gamma]] := Module[
  { n = Length[p], foc, leqs},
  (* Get the derivatives *)
  foc = Diagonal[dvector[p, cost, v, g]];
  (* Set the equations *)
  leqs = Thread[foc == ConstantArray[0, n]];
  (* Solve the equations *)
  Solve[leqs, p][[-1]]
  ]

(* Check on an example *)
solution[{x, y, z}]

Now, the solution function is fully functionnal on any example price vector of defined length like {x,y,z} or {p1, p2, p3, p4, p5} for instance.

Yet, my goal is to find the general solution of the system, which is a function of N the length of my price vector. More precisely, if p(1),...,p(N) is my price vector, I would like to express, for any i, p(i) as a function of the primitives of the model as well as N, the length of the price vector.

Is there a way to solve this in its general form ?

Thanks in advance,

M

Answer 1

In Mathematica scalars eg $v$ and vectors can be mixed gracefully, courtesy of the beautiful kernel; this means, practically, that if eg $p$ is a vector or-more loosely speaking (this, too, is a distinction that is sometimes blurred but certainly not without utility)-a List, then something like v+p evaluates into the equivalent of $(v+p_1,...,v+p_n)$.

Also, when summing over a list, Total is a more suitable choice instead of Sum; the later makes use of indexes (the i part) which is not needed in order to access a list. Mathematica, unlike traditional procedural languages, does not need (it surely can, but it's not compelled) to use indexing of lists. In the present context, it wouldn't do much harm I guess but having to deal with the extra indexes makes it a bad practice to sum over immutable lists like that.

So, having said that, the following is an implementation of quantite:

  Clear[V, γ]
  q[p_, v_: V, g_: γ] := Module[{meanP = Mean[p]},
    V - p - \[Gamma] ( p - meanP )
   ]

The argument of q is intended to be a vector of prices and the output is a vector of quantities. Evaluate q[{a,b}] to verify that for the price vector {a, b} the output is as expected.

Similarly, define the profit function as

Clear[c];
profit[p_, const_:c, v_: V, g_: γ] := (p - const) q[p, v, g]

Again, evaluate profit[{a, b}] to make sure the output is appropriate for input price vector {a, b}.

Finally, formulating the foc equations does not need to make use of indexing.

dvector[p_, const_: c, v_: V, g_: γ] := D[profit[p, const, v, g], {p}]

Admittedly, this implementation of dvector is not the most efficient; the derivatives (actually, it's the Jacobian of profit) are recalculated for every different price vector; a more robust implementation would probably first calculate the derivatives and store the result for later use; the chosen implementation is coded faster and is more transparent to someone with little experience with Mathematica.

At this point, please make sure you understand how D does indeed produce the Jacobian.

Solving for the prices is fairly straightforward:

solution[p_, const_: c, v_: V, g_: γ] := Module[{n = Length[p], foc, veqs},
  (* obtain the derivatives *)
  foc = Diagonal[dvector[p, const, v, g]];
  (* set-up the equations *)
  veqs = Thread[foc == ConstantArray[0, n]];

  Solve[veqs, p][[-1]]
 ]

Make sure to evaluate solution[{a, b}] to verify the solution is as expected.