Skip to content

Commit 2937564

Browse files
SebKrantzSebKrantz
committed
First version that is almost working.
1 parent 482f415 commit 2937564

1 file changed

Lines changed: 63 additions & 104 deletions

File tree

‎src/models/solve_allocation_by_duality_cgc.jl‎

Lines changed: 63 additions & 104 deletions
Original file line numberDiff line numberDiff line change
@@ -142,13 +142,19 @@ function hessian_duality_cgc(
142142
Zjn = graph.Zjn
143143
Lj = graph.Lj
144144
omegaj = graph.omegaj
145-
# N = param.N
146145
hess_str = auxdata.hess
147146

147+
# Constant term: first part of Bprime
148+
cons = m1dbeta * n1dnum1 * numbetam1
149+
# Constant in T'
150+
cons2 = numbetam1 / (numbetam1+nu*beta)
151+
148152
# Precompute elements
149153
res = recover_allocation_duality_cgc(x, auxdata)
150154
Pjn = res.Pjn
151-
PDj = res.PDj
155+
PCj = res.PCj
156+
Qjk = res.Qjk
157+
Qjkn = res.Qjkn
152158

153159
# Prepare computation of production function
154160
if a != 1
@@ -173,104 +179,57 @@ function hessian_duality_cgc(
173179
# Stores the derivative term
174180
term = 0.0
175181

176-
# Needed for both Q^n_{jk} and D^n_j
177-
P = PDj[j]^m1dbeta
178-
Pprime = m1dbeta * Pjn[j, nd]^(-sigma) * P^(1+beta - beta*sigma)
179-
180-
# Starting with D^n_j = (C+T)G, derivative: (C'+T')G + (C+T)G'
181-
# T' can be additively split. Here I compute (C' + T'_1)G + (C+T)G'
182-
G = (Pjn[j, n] / PDj[j])^(-sigma)
182+
# Starting with D^n_j = (C+T)G, derivative: (C'+T')G + (C+T)G' = C'G + CG' + T'G + TG'
183+
# Terms involving T are computed below. Here forcus on C'G + CG'
184+
G = (Pjn[j, n] / PCj[j])^(-sigma)
183185
if jd == j # 0 for P^n_k
184-
T = 0.0
185-
for k in neighbors # TODO: add to loop below ??
186-
T += res.Qjk[j, k]^(1+beta) / kappa[j, k]
187-
end
188-
Gprime = sigma * (Pjn[j, n] * Pjn[j, nd])^(-sigma) * PDj[j]^(2*sigma-1)
189-
# Adding all derivative terms apart from second part of T'
186+
Gprime = sigma * (Pjn[j, n] * Pjn[j, nd])^(-sigma) * PCj[j]^(2*sigma-1)
187+
Cprime = Lj[j]/omegaj[j] * (PCj[j] / Pjn[j, nd])^sigma / param.usecond(res.cj[j], graph.hj[j]) # * param.uprimeinvprime(PCj[j]/omegaj[j], graph.hj[j])
190188
if nd == n
191-
Cprime = Lj[j]/omegaj[j] * (PDj[j] / Pjn[j, nd])^sigma * param.uprimeinvprime(PDj[j]/omegaj[j], graph.hj[j])
192-
term += (Cprime + (1+beta) * Pprime / P * T) * G
193-
Gprime -= sigma / PDj[j] * G^((sigma+1)/sigma)
189+
Gprime -= sigma / PCj[j] * G^((sigma+1)/sigma)
194190
end
195-
term += (res.Cj[j] + T) * Gprime
191+
term += Cprime * G + res.Cj[j] * Gprime
196192
end
197-
# Constant term for iterating through neighbors for computing T'_2 G
198-
cons = numbetam1 / (numbetam1 + nu*beta) * P^(1+beta) * G
199-
# Cconstant tern: first part of Bprime
200-
cons2 = m1dbeta * n1dnum1 * numbetam1
201193

202-
# Now terms sum_k(Q^n_{jk} - Q^n_{kj}) as well as T'_2 G
194+
# Now terms sum_k(Q^n_{jk} - Q^n_{kj}) as well as T'G + TG'
203195
for k in neighbors
204-
205-
# Terms with Derivative
206-
diff = Pjn[k, nd] - Pjn[j, nd]
207-
if diff >= 0
208-
A = (diff/m[nd])^(1/(nu-1)) # 0 for n' != n
209-
# Constant Terms
210-
K0 = (1 + beta) / kappa[j, k]
211-
PK0 = PDj[j] * K0
212-
K = K0^m1dbeta
213-
214-
# Terms with Derivatives..continued
215-
B = (res.Qjk[j, k] * PK0^m1dbeta)^(numbetam1/(nu-1))
216-
217-
# Computing the right derivative: Q^n_{jk}
218-
if jd == j # P^x_j
219-
Bprime = cons2 * A * B^((numbetam1 - nu*beta)/numbetam1)
220-
if nd == n # P^n_j
221-
Aprime = -n1dnum1 / m[n] * A^(2-nu)
222-
term += K * (Pprime*A*B + P*Aprime*B + P*A*Bprime)
223-
else # P^n'_j (A' is zero)
224-
term += K * (Pprime*A*B + P*A*Bprime)
225-
end
226-
# This computes the remaining derivative of D^n_j (T'_2 G)
227-
term += cons * K / (1+beta) * Bprime * B^((-nu*beta)/(numbetam1+nu*beta))
228-
elseif jd == k # P^x_k (P' is zero)
229-
Bprime = -cons2 * A * B^((numbetam1 - nu*beta)/numbetam1) # simply the negative
230-
if nd == n # P^n_k
231-
Aprime = n1dnum1 / m[n] * A^(2-nu) # Simply the negative
232-
term += K * (P*Aprime*B + P*A*Bprime)
233-
else # # P^n'_k (A' is zero)
234-
term += K * P*A*Bprime
235-
end
236-
# This computes the remaining derivative of D^n_j (T'_2 G)
237-
term += cons * K / (1+beta) * Bprime * B^((-nu*beta)/(numbetam1+nu*beta))
196+
diff = Pjn[k, n] - Pjn[j, n]
197+
if diff >= 0 # Flows in the direction of k
198+
tmp = Qjkn[j, k, n]
199+
PK0 = PCj[j] * (1 + beta) / kappa[j, k]
200+
KPABprimemQjkn = cons * ((Pjn[k, nd] - Pjn[j, nd])/m[nd])^n1dnum1 * Qjk[j, k]^(-nu*beta*n1dnum1) * PK0^(nu*n1dnum1)
201+
if jd == j
202+
KPprimeABmQjkn = m1dbeta * Pjn[j, nd]^(-sigma) * PCj[j]^(sigma-1)
203+
term += tmp * KPprimeABmQjkn # KP'AB
204+
if nd == n
205+
term -= tmp * n1dnum1 / diff # KPA'B
206+
end
207+
term += tmp * KPABprimemQjkn # KPAB'
208+
term += Qjk[j, k]^(1+beta) / kappa[j, k] * ((1+beta)*KPprimeABmQjkn + cons2 * KPABprimemQjkn) * G # T'G
209+
term += Qjk[j, k]^(1+beta) / kappa[j, k] * Gprime # TG'
210+
elseif jd == k
211+
if nd == n
212+
term += tmp * n1dnum1 / diff # KPA'B [A'(k) has opposite sign]
213+
end
214+
term -= tmp * KPABprimemQjkn # KPAB' [B'(k) has opposite sign]
215+
term -= Qjk[j, k]^(1+beta) / kappa[j, k] * cons2 * KPABprimemQjkn * G # T'G [B'(k) has opposite sign]
238216
end
239-
end
240-
241-
# Computing the right derivative: Q^n_{kj} (other direction -> needs to be subtracted)
242-
243-
# Terms with Derivatives
244-
if diff < 0
245-
A = (-diff/m[nd])^(1/(nu-1)) # 0 for n' != n
246-
# Also need prices here
247-
Pk = PDj[k]^m1dbeta
248-
Pkprime = m1dbeta * Pjn[k, nd]^(-sigma) * Pk^(1+beta - beta*sigma)
249-
250-
# Constant Terms
251-
K0 = (1 + beta) / kappa[k, j]
252-
PK0 = PDj[k] * K0
253-
K = K0^m1dbeta
254-
255-
# Terms with Derivatives..continued
256-
B = (res.Qjk[k, j] * PK0^m1dbeta)^(numbetam1/(nu-1))
257-
# Computing the right derivative: Q^n_{kj}
258-
if jd == k # P^x_k
259-
Bprime = cons2 * A * B^((numbetam1 - nu*beta)/numbetam1)
260-
if nd == n # P^n_k
261-
Aprime = -n1dnum1 / m[n] * A^(2-nu)
262-
term -= K * (Pkprime*A*B + Pk*Aprime*B + Pk*A*Bprime)
263-
else # P^n'_k (A' is zero)
264-
term -= K * (Pkprime*A*B + Pk*A*Bprime)
265-
end
266-
elseif jd == j # P^x_j (P' is zero)
267-
Bprime = -cons2 * A * B^((numbetam1 - nu*beta)/numbetam1) # simply the negative
268-
if nd == n # P^n_j
269-
Aprime = n1dnum1 / m[n] * A^(2-nu) # Simply the negative
270-
term -= K * (Pk*Aprime*B + Pk*A*Bprime)
271-
else # # P^n'_j (A' is zero)
272-
term -= K * Pk*A*Bprime
273-
end
217+
else # Flows in the direction of j
218+
tmp = Qjkn[k, j, n]
219+
PK0 = PCj[k] * (1 + beta) / kappa[k, j]
220+
KPABprimemQjkn = cons * ((Pjn[j, nd] - Pjn[k, nd])/m[nd])^n1dnum1 * Qjk[k, j]^(-nu*beta*n1dnum1) * PK0^(nu*n1dnum1)
221+
if jd == k
222+
KPprimeABmQjkn = m1dbeta * Pjn[k, nd]^(-sigma) * PCj[k]^(sigma-1)
223+
term -= tmp * KPprimeABmQjkn # KP'AB
224+
if nd == n
225+
term += tmp * n1dnum1 / abs(diff) # KPA'B
226+
end
227+
term -= tmp * KPABprimemQjkn # KPAB'
228+
elseif jd == j
229+
if nd == n
230+
term -= tmp * n1dnum1 / abs(diff) # KPA'B [Akj'(j) has opposite sign]
231+
end
232+
term += tmp * KPABprimemQjkn # KPAB' [Bkj'(j) has opposite sign]
274233
end
275234
end
276235
end # End of k loop
@@ -285,7 +244,7 @@ function hessian_duality_cgc(
285244
end
286245

287246
# Assign result
288-
values[ind] = -obj_factor * term
247+
values[ind] = -obj_factor * term + 1e-5
289248
end
290249
end
291250
return
@@ -305,7 +264,7 @@ function recover_allocation_duality_cgc(x, auxdata)
305264

306265
# Extract price vectors
307266
Pjn = reshape(x, (graph.J, param.N))
308-
PDj = sum(Pjn .^ (1 - param.sigma), dims=2) .^ (1 / (1 - param.sigma))
267+
PCj = sum(Pjn .^ (1 - param.sigma), dims=2) .^ (1 / (1 - param.sigma))
309268

310269
# Calculate labor allocation
311270
if param.a < 1
@@ -320,14 +279,14 @@ function recover_allocation_duality_cgc(x, auxdata)
320279
Yjn = graph.Zjn .* (Ljn .^ param.a)
321280

322281
# Calculate aggregate consumption
323-
cj = param.alpha * (PDj ./ omegaj) .^
282+
cj = param.alpha * (PCj ./ omegaj) .^
324283
(-1 / (1 + param.alpha * (param.rho - 1))) .*
325284
hj1malpha .^ (-(param.rho - 1)/(1 + param.alpha * (param.rho - 1)))
326285
Cj = cj .* Lj
327286

328287
# Calculate the aggregate flows Qjk
329288
Qjk = zeros(graph.J, graph.J)
330-
temp = (kappa ./ ((1 + param.beta) * PDj)) .^ (nu/(nu-1))
289+
temp = (kappa ./ ((1 + param.beta) * PCj)) .^ (nu/(nu-1))
331290
temp[nadj] .= 0
332291
for n in 1:param.N
333292
Lambda = repeat(Pjn[:, n], 1, graph.J)
@@ -338,7 +297,7 @@ function recover_allocation_duality_cgc(x, auxdata)
338297

339298
# Calculate the flows Qjkn
340299
Qjkn = zeros(graph.J, graph.J, param.N)
341-
temp = kappa ./ ((1 + param.beta) * PDj .* Qjk .^ (1+param.beta-nu))
300+
temp = kappa ./ ((1 + param.beta) * PCj .* Qjk .^ (1+param.beta-nu))
342301
temp[nadj .| isinf.(temp)] .= 0 # Because of the max() clause, Qjk may be zero
343302
for n in 1:param.N
344303
Lambda = repeat(Pjn[:, n], 1, graph.J)
@@ -349,9 +308,9 @@ function recover_allocation_duality_cgc(x, auxdata)
349308
temp = Qjk .^ (1+param.beta) ./ kappa
350309
temp[nadj] .= 0
351310
Dj = Cj + sum(temp, dims = 2)
352-
Djn = Dj .* (Pjn ./ PDj) .^ (-param.sigma)
311+
Djn = Dj .* (Pjn ./ PCj) .^ (-param.sigma)
353312

354-
return (Pjn=Pjn, PDj=PDj, Ljn=Ljn, Yjn=Yjn, cj=cj, Cj=Cj, Dj=Dj, Djn=Djn, Qjk=Qjk, Qjkn=Qjkn)
313+
return (Pjn=Pjn, PCj=PCj, Ljn=Ljn, Yjn=Yjn, cj=cj, Cj=Cj, Dj=Dj, Djn=Djn, Qjk=Qjk, Qjkn=Qjkn)
355314
end
356315

357316

@@ -382,10 +341,10 @@ end
382341
# mat_kappa = repeat(kappa, param.N, param.N)
383342
# mN = repeat(param.m, inner = graph.J)
384343
# # Rows are the derivatives, columns are P^n_j
385-
# PDjN = repeat(res.PDj, param.N)
386-
# A = mat_kappa ./ ((1+param.beta) * mN .* PDjN)'
344+
# PCjN = repeat(res.PCj, param.N)
345+
# A = mat_kappa ./ ((1+param.beta) * mN .* PCjN)'
387346
# # Adding term for (block-digonal) elements
388-
# temp = diff .* (x .^ (-param.sigma) .* PDjN .^ (param.sigma - 1))' .- 1
347+
# temp = diff .* (x .^ (-param.sigma) .* PCjN .^ (param.sigma - 1))' .- 1
389348
# temp[Iij .== 0] .= 1
390349
# Aprime = A .* temp
391350
# temp = diff .* Inm
@@ -406,7 +365,7 @@ end
406365
# # Derivatives of Qjk
407366
# # Result should be nedg * N
408367
# Qjkprime = zeros(graph.J, graph.J, param.N)
409-
# temp = kappa ./ ((1 + param.beta) * res.PDj)
368+
# temp = kappa ./ ((1 + param.beta) * res.PCj)
410369
# for n in 1:param.N
411370
# Lambda = repeat(Pjn[:, n], 1, graph.J)
412371
# Qjkprime[:,:,n] = m[n] / param.beta * res.Qjk .^ ((nu-nu*param.beta-1)/(nu-1))

0 commit comments

Comments
 (0)