Hello I have this modified monocentric city model that I have been working on for a while.
Sorry for the long post, I wanted to add the whole code, because it's all connected.
My suspicion is that it has something to do with the p (price) and q (demand) functions [economics]. My previous version used a simpler price and demand functions, but since I changed them, I have not been able to get it to work. I am no mathematical genius, but I have tried so many things that I'm about to give up hope on this model. I hope one of you can help me figure out the problem.
(* Cost share of housing capital in housing production function *)
\[Beta] = 0.6;
(* Scaling on housing production function *)
g = 0.0005;
(* Radians available for construction *)
(* benchmark \[Theta] is 3 *)
\[Theta] = 2*Pi;
(* elasticity of income *)
\[Gamma] = 0.2;
(* price elasticity of housing *)
\[Epsilon] = -0.6;
(* constant in demand function *)
\[Omega] = 1;
p[x_, y_, t_, ta_, f_, fa_,
u_] := (E^\[Omega]*(\[Epsilon] +
1)*((y - (f + fa) - (t + ta)*x)^(1 - \[Gamma])/(1 - \[Gamma]) -
u))^(1/(\[Epsilon] + 1));
q[x_, y_, t_, ta_, f_, fa_,
u_] := (1 - \[Gamma])^(\[Gamma]/(1 - \[Gamma]))*E^\[Omega]*
p[x, y, t, ta, f, fa,
u]^(\[Epsilon] +
1)*((E^\[Omega]*
p[x, y, t, ta, f, fa, u]^(\[Epsilon] + 1))/(\[Epsilon] + 1) +
u)^(\[Gamma]/(1 - \[Gamma]));
S[x_, y_, t_, ta_, f_, fa_,
u_] := (1/(p[x, y, t, ta, f, fa, u] (\[Beta]) (g)))^(1/(\[Beta] -
1));
r[x_, y_, t_, ta_, f_, fa_,
u_] := (p[x, y, t, ta, f, fa,
u] (g)) ((1/(p[x, y, t, ta, f, fa,
u] (\[Beta]) (g)))^(\[Beta]/(\[Beta] - 1))) -
1 ((1/(p[x, y, t, ta, f, fa, u] (\[Beta]) (g)))^(1/(\[Beta] - 1)));
h[x_, y_, t_, ta_, f_, fa_,
u_] := (g) S[x, y, t, ta, f, fa, u]^(\[Beta]);
Density[x_, y_, t_, ta_, f_, fa_, u_] :=
h[x, y, t, ta, f, fa, u]/q[x, y, t, ta, f, fa, u];
L[xbar_, y_, t_, ta_, f_, fa_, u_] := \[Theta]*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*\
Density[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
pbar[xbar_, y_, t_, ta_, f_, fa_, u_,
pop_] := (\[Theta]/pop)*(1/30000)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*
Density[x, y, t, ta, f, fa, u]*
p[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
qbar[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*
Density[x, y, t, ta, f, fa, u]*
q[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
hbar[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x*
Density[x, y, t, ta, f, fa, u]*
h[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
xhat[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] := (\[Theta]/pop)*\!\(
\*SubsuperscriptBox[\(\[Integral]\), \(0\), \(xbar\)]\(x^2*
Density[x, y, t, ta, f, fa, u] \[DifferentialD]x\)\);
rev[xbar_, y_, t_, ta_, f_, fa_, u_, pop_] :=
pop*(fa + ta*xhat[xbar, y, t, ta, f, fa, u, pop]);
(*fa[xbar_,y_,t_,ta_,f_,fa_,u_,pop_] := \
((rev[xbar,y,t,ta,f,fa,u,pop]/pop)/(ta*xhat[xbar,y,t,ta,f,fa,u,pop]));
xhatnew[xbar_,y_,t_,ta_,f_,fa_,u_,pop_]:= \
((rev[xbar,y,t,ta,f,fa,u,pop]/pop)-fa[xbar,y,t,ta,f,fa,u,pop])/ta;*)
Population = 800000;
y = 70000;
t = 350;
f = 3000;
ta = 250;
fa = 1000;
ra = 45000;
Lu[u_] := (solution =
FindRoot[{r[xbar, y, t, ta, f, fa, u] - ra == 0}, {xbar, 10},
AccuracyGoal -> 1]; xbar /. solution);
fu[u_] := N[L[Evaluate[Lu[u]], y, t, ta, f, fa, u] - Population];
findRoot[fun_, lower_, upper_, eps_] :=
Module[{ul = lower, uh = upper, count = 0},
While[And[Abs[fun[ul]] > eps, count < 10],
If[fun[(uh + ul)/2] > 0, ul = (uh + ul)/2, uh = (uh + ul)/2];
count = count + 1]; Print["Iterations: ", count,
" Population: ", Population, " people",
" utility: ", N[ul],
" City boundary (xbar): ", Lu[N[ul]], " km",
" Average house price (pbar): ",
pbar[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " euros",
" Price at 5 km from the centre p(5): ",
p[5, y, t, ta, f, fa, N[ul]]/30000, " euros",
" Price at 10 km from the centre p(10): ",
p[10, y, t, ta, f, fa, N[ul]]/30000, " euros",
" Average density: ", 2*Population/(\[Theta]*Lu[N[ul]]^2),
" people per ?",
" Average floor to area ratio hbar: ",
hbar[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population],
" Density at 5 km from the centre Density(5): ",
Density[5, y, t, ta, f, fa, N[ul]], " people per ?",
" Density at 20 km from the centre Density(20): ",
Density[20, y, t, ta, f, fa, N[ul]], " people per ?",
" Average distance to the centre xhat: ",
xhat[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " km",
" Tax revenue from transport: ",
rev[Lu[N[ul]], y, t, ta, f, fa, N[ul], Population], " euros",
Plot[{Density[x, y, t, ta, f, fa, N[ul]]}, {x, 0, Lu[N[ul]]}]]];
findRoot[fu, 3000, 10000, 1];
It just keeps running. I don't know how long, cause after a day I stopped it.
Extra details:
The p is derived from this indirect utility function (log utility as opposed to my previous Cobb-Douglas utility) through simple equation solving assuming that u is constant:
u = (-e^(\[Omega]*∙))*(p^(\[Epsilon] + 1)/(\[Epsilon] + 1)) +
y^(-\[Gamma] + 1)/(-\[Gamma] + 1)
y is to be replaced by y-(fa+f)-(ta+t)x
The q is derived by taking the u functionand putting y-(fa+f)-(ta+t)x on the left side:
y - (f + fa) - (t + ta)*
x = ((-\[Gamma] + 1)*(u +
e^(\[Omega]*∙)*(p[x, y, t, ta, f, fa,
u]^(\[Epsilon] + 1)/(\[Epsilon] + 1))))^(1/(-\[Gamma] + 1))
and differentiating with respect to p, which should be the equation above for q.
Any help would be greatly appreciated. Also help about the process of troubleshooting this.
