@@ -320,97 +320,3 @@ function recover_allocation_duality_cgc(x, auxdata)
320320 return (Pjn= Pjn, PCj= PCj, Ljn= Ljn, Yjn= Yjn, cj= cj, Cj= Cj, Dj= Dj, Djn= Djn, Qjk= Qjk, Qjkn= Qjkn)
321321end
322322
323-
324-
325- # # Hessian Experimental:
326- # # These are the lagrange multipliers = prices
327- # Lambda = repeat(x, 1, graph.J * param.N) # P^n_j: each column is a price, the rows should be the derivatives (P^n'_k)
328- # lambda = reshape(x, (graph.J, param.N))
329-
330- # # Compute price index
331- # P = sum(lambda .^ (1 - param.sigma), dims=2) .^ (1 / (1 - param.sigma))
332- # mat_P = repeat(P, param.N, param.N * graph.J)
333-
334- # # Create masks
335- # # This lets different products in the same location relate
336- # Iij = kron(ones(param.N, param.N), I(graph.J))
337- # # This is adjacency, but only for the same product (n == n')
338- # Inm = kron(I(param.N), graph.adjacency)
339-
340- # # Compute Qjkn terms for Hessian
341- # Qjknprime = zeros(graph.J, graph.J, param.N)
342-
343-
344-
345-
346-
347- # diff = Lambda' - Lambda # P^n'_k - P^n_j
348- # mat_kappa = repeat(kappa, param.N, param.N)
349- # mN = repeat(param.m, inner = graph.J)
350- # # Rows are the derivatives, columns are P^n_j
351- # PCjN = repeat(res.PCj, param.N)
352- # A = mat_kappa ./ ((1+param.beta) * mN .* PCjN)'
353- # # Adding term for (block-digonal) elements
354- # temp = diff .* (x .^ (-param.sigma) .* PCjN .^ (param.sigma - 1))' .- 1
355- # temp[Iij .== 0] .= 1
356- # Aprime = A .* temp
357- # temp = diff .* Inm
358- # temp[Inm .== 0] .= 1
359- # A .*= temp
360-
361- # BA = A .^ (1/(nu-1))
362- # BprimeAAprime = (1/(nu-1)) * A .^ (nu/(1-nu)) .* Aprime
363- # BprimeAAprime[isnan.(BprimeAAprime)] .= 0
364-
365- # C = repeat(res.Qjk .^ ((nu-param.beta-1)/(nu-1)), param.N, param.N)
366- # Cprime = repeat(((nu-param.beta-1)/(nu-1)) * res.Qjk .^ (param.beta/(1-nu)), param.N, param.N)
367-
368- # Qjknprime = BA .* Cprime .* BprimeAAprime .* mN # Off-diagonal (n') part
369- # Qjknprime[Inm] += BprimeAAprime[Inm] .* C[Inm] # Adding diagonal
370- # # Compute Qjkn Sums
371-
372- # # Derivatives of Qjk
373- # # Result should be nedg * N
374- # Qjkprime = zeros(graph.J, graph.J, param.N)
375- # temp = kappa ./ ((1 + param.beta) * res.PCj)
376- # for n in 1:param.N
377- # Lambda = repeat(Pjn[:, n], 1, graph.J)
378- # Qjkprime[:,:,n] = m[n] / param.beta * res.Qjk .^ ((nu-nu*param.beta-1)/(nu-1))
379- # LL = Lambda' - Lambda # P^n_k - P^n_j
380- # LL[.!graph.adjacency] .= 0
381- # Qjkprime[:,:,n] .*= ((LL .* temp) / m[n]) .^ (1/(nu-1)) # A^(1/(nu-1))
382- # Qjkprime[:,:,n] .*= temp / m[n] # A'(P^n_k)
383- # tril()
384- # end
385- # Qjk = (Qjk .* temp) .^ (nu-1)/(nu*param.beta)
386-
387- # Qjkprime = res.Qjk[graph.adjacency]
388-
389-
390- # # This is for Djn, the standard trade part (excluding trade costs)
391- # termA = -param.sigma * (repeat(P, param.N) .^ param.sigma .* x .^ (-(param.sigma + 1)) .* repeat(res.Cj, param.N))
392- # part1 = Iij .* Lambda .^ (-param.sigma) .* Lambda' .^ (-param.sigma) .* mat_P .^ (2 * param.sigma)
393- # termB = param.sigma * part1 ./ mat_P .* repeat(res.Cj, param.N, graph.J * param.N)
394- # termC = part1 .* repeat(graph.Lj ./ (graph.omegaj .* param.usecond.(res.cj, graph.hj)), param.N, graph.J * param.N)
395- # Cjn = Diagonal(termA[:]) + termB + termC
396-
397- # # Now comes the part of Djn relating to trade costs
398- # Djn_costs =
399-
400-
401- # Djn = Cjn + Djn_costs
402-
403-
404- # # This is related to the flows terms in the standard case
405-
406- # # Off-diagonal P^n_k terms
407- # termD = 1 / (param.beta * (1 + param.beta)^(1 / param.beta)) * Inm .* mat_kappa .^ (1 / param.beta) .*
408- # abs.(diff) .^ (1 / param.beta - 1) .*
409- # ((diff .> 0) .* Lambda' ./ Lambda .^ (1 + 1 / param.beta) +
410- # (diff .< 0) .* Lambda ./ Lambda' .^ (1 + 1 / param.beta))
411- # # Diagonal P^n_j terms: sum across kappa
412- # termE = -1 / (param.beta * (1 + param.beta)^(1 / param.beta)) * Inm .* mat_kappa .^ (1 / param.beta) .*
413- # abs.(diff) .^ (1 / param.beta - 1) .*
414- # ((diff .> 0) .* Lambda' .^ 2 ./ Lambda .^ (2 + 1 / param.beta) +
415- # (diff .< 0) ./ Lambda' .^ (1 / param.beta))
416- # termE = sum(termE, dims=2)
0 commit comments