$title Pilot Implementation of Bernard, Redding and Schott (Version 5)

$ontext

Reference: "Comparative Advantage and Heterogeneous Firms", Working
Paper, July 2004 by Andrew B. Bernard, Stephen Redding and Peter K.
Schott.

I found the model formulated in this paper to be very interesting, but
computationally infeasible. The problem is that the definition of an
equilibrium in a model with N regions and M markets involves NxNxM
entry decisions for monopolistic firms. The curse of dimensionality
suggests that while the results are intriguing, it may be impossible
to evaluate the model in an empirical context.


The file illustrates a decomposition method which provides some hope
for empirical applications of this framework.

I propose a decomposition through which we can solve multiregional,
multisectoral general equilibrium models with heterogeneous firms
engaging in monopolistic competitive markets. A partial equilibrium
model of firm-level, multi-regional monopolistic competition is
embedded within a general equilibrium model. The general equilibrium
model determines factor prices and comparative advantage. The partial
equilibrium model assesses productivity effects which arise from the
introduction of new varieties within a Dixit-Stiglitz demand
framework.

In the partial equilibrium model firm entry and exit decisions are
made on the basis of expected profitability, depending on productivity
levels drawn from a Pareto distribution. Exporting activities are
determined by firm-specific productivity levels. Only those firms
with sufficiently high productivity engage in trade, as trade flows
are subject to fixed costs.

The model incorporates Dixit-Stiglitz preferences in final demand with
firm-level product differentiation. Factor markets are competitive
and introduce conventional comparative advantage factors into the
model.

This version includes an upward sloping supply function for firm-level
resources in the partial equilibrium subproblem. The slope of the
supply function is taken as 5, although this parameter need not be
exact.

$offtext

set     i Industries /1,2/
        k Regions /H,F/
        f Primary factors /
                s Skilled labor,
                u Unskilled labor/,
        t(i,k,k) Set of active trade links,
        eq(i) Active industries;

eq(i) = yes;

alias (k,kk,kkk);

parameter

* Inputs describing Pareto distribution for firm productivity

*       g(phi) = (b/phi)**a * a/phi (36)
* 
*       a > (sigma-1) so that log firm sales have a finite variance

* Standard deviation of the log of firm sales is 1/(a-sigma+1)

        a       Shape parameter in the Pareto distribution /3.4/
        b       Minimum productivity in the Pareto distribution /0.2/
        sigma   Intra-variety elasticity /3.8/
        beta(i,f) Factor intensities /1.s 0.6, 2.s 0.4/
        alpha(i) Value shares /1 0.5, 2 0.5/
        fe(i)   Fixed cost of entry /1 2, 2 2/
        fx      Fixed costs of export /0.2/
        delta   Probability of firm death /0.025/,
        eta     Elasticity of supply in PE submodel /5/
        rho     Inverse elasticity
        fc(i,k,k) Bilateral fixed cost,
        tau(i,k,k) Iceberg trade cost parameter;

fc(i,k,kk) = fx;
fc(i,k,k) = 0.1;

*       Begin with autarchy and then explore tau ranges from 2 to 1.2

*       tau(i,k,kk) = +inf;

tau(i,k,kk) = 2;
tau(i,k,k) = 1;
t(i,k,kk) = yes$(1/tau(i,k,kk));
beta(i,"u") = 1 - beta(i,"s");
rho = 1 - 1/sigma;

*       Parameters calibrated to match plant-level manufacturing data
*       in Bernard et al (2003)

table s(f,k) Primary factor endowments
        H       F
S       1200    1000
U       1000    1200;

VARIABLES
        P(i,k)          Price index
        PT(i,k,kk)      Weighted average price index
        QT(i,k,kk)      Mean bilateral sales
        M(i,k)          Total number of firms
        N(i,k,kk)       Number of firms
        THETA(i,k,kk)   Firm supply shares
        PHIS(i,k,kk)    Minimum productivity,
        PI(i,k,kk)      Bilateral profit
        PHIT(i,k,kk)    Weighted average productivity,
        Y(i,k)          Aggregate sectoral output,
        C(i,k)          Unit cost of firm inputs
        W(f,k)          Factor price
        R(k)            Regional income;

POSITIVE VARIABLES N, M, QT, W, P, C, Y, R;

EQUATIONS ptdef, pdef, qtdef, phisdef, pidef, phitdef, thetadef,
ndef, mdef, ydef, cdef, cdefpe, fmarket, income;

*       The unit cost of goods are defined on the basis of cost-minimizing
*       behavior with Cobb-Douglas technology:
cdef(i,k).. C(i,k) =e= prod(f, W(f,k)**beta(i,f));

*       Factor supply equals facgtor demand (the market clearance equations
*       are multiplied by the factor price to improve numerical stability):

fmarket(f,k).. W(f,k)*s(f,k) =e= sum(i, beta(i,f)*C(i,k)*Y(i,k));

*       Factor income is based on factor prices and endowments:

income(k).. R(k) =e= sum(f, W(f,k)*s(f,k));

*       Aggregate sectoral output is based on the number of firms. Sectoral
*       output consists of fixed costs of supply, fixed trade costs and
*       variable costs of suply which depend on the quantity shipped, the
*       trade costs and weighted average productity:

ydef(i,k)$eq(i)..
        Y(i,k) =e= M(i,k)* (delta*fe(i) + sum(t(i,k,kk),
                THETA(i,k,kk) * (fc(i,k,kk)+QT(i,k,kk)*tau(i,k,kk)/PHIT(i,k,kk))));

*       In the partial equilibrium model we fix industry output but assume that
*       supply responds to shadow value of output. This equation has an
*       exception operator eq(i) indicating that it only enters for the active
*       sector:

cdefpe(i,k)$eq(i)..
        Y.L(i,k)*(C(i,k)/C.L(i,k))**eta =e= M(i,k)* (delta*fe(i) + 
          sum(t(i,k,kk), THETA(i,k,kk) * (fc(i,k,kk)+QT(i,k,kk)*tau(i,k,kk)/PHIT(i,k,kk))));

*       Free entry zero profit condition relates fixed cost of entry
*       to profits and -- this determines the aggregate number of firms
*       operating in region k:

mdef(i,k)$eq(i)..  C(i,k)*fe(i) =e= sum(t(i,k,kk), PI(i,k,kk) * THETA(i,k,kk))/delta;

*       Output per firm is constant in this model (due to the
*       assumption of iceberg trade cost and identical factor content
*       of fixed and variable costs), so the supply from kk to k
*       is the same as the number of firms:

pdef(i,k)$eq(i)..  sum(t(i,kk,k), N(i,kk,k)*PT(i,kk,k)**(1-sigma)) =e= P(i,k)**(1-sigma);

*       The cost of supply from k to kk depends on factor cost (C), the
*       transportation markup (tau), the weighted average productivity
*       (PHIT) and the markup factor (rho):

ptdef(t(i,k,kk))..      PT(i,k,kk) * (rho * PHIT(i,k,kk)) =e= tau(i,k,kk) * C(i,k);

*       The aggregate sales from k into kk follows from a nested
*       CES demand function in which the top-level preferences are
*       are Cobb-Douglas:
qtdef(t(i,k,kk))..     QT(i,k,kk) =e= (alpha(i)*R(kk)/P(i,kk)) * (P(i,kk)/PT(i,k,kk))**sigma;

*       Unit profits on individual firm sales from k to kk equal
*       markup revenue less fixed costs:

pidef(t(i,k,kk))..      PI(i,k,kk) =e= PT(i,k,kk)*QT(i,k,kk)/sigma - C(i,k)*fc(i,k,kk);

* The fraction of firms from market k selling into market kk:

thetadef(t(i,k,kk))..   THETA(i,k,kk) =e= (b/PHIS(i,k,kk))**a;

*       The number of firms selling from k to kk depend on the number of 
*       firms operating in region k:

ndef(t(i,k,kk))..       N(i,k,kk) =e= THETA(i,k,kk) * M(i,k);

*       The minimum productivity level firm selling from region k
*       to region kk:

phisdef(t(i,k,kk))..    QT(i,k,kk)*PT(i,k,kk) * ((a+1-sigma)/a) =e= sigma*fc(i,k,kk)*C(i,k);

*       The weighted average productivity is computed on the basis of
*       the cutoff productivity level for firms in the industry.

phitdef(t(i,k,kk))..    PHIT(i,k,kk) =e= (a/(a+1-sigma))**(1/(sigma-1)) * PHIS(i,k,kk);

*       Define two models, one which includes the full systems and a second
*       which defines a partial equilibrium model for a single commodity:

MODEL   MCHF /  pdef.P, ndef.N, ptdef.PT, qtdef.QT, pidef.PI, thetadef.THETA,
                phisdef.PHIS, phitdef.PHIT, ydef.Y,
                mdef.M, cdef.C, fmarket.W, income.R /,

        MCHFPE /pdef.P, ndef.N, ptdef.PT, qtdef.QT, pidef.PI, thetadef.THETA,
                phisdef.PHIS, phitdef.phit, mdef.M, cdefpe.C/;

*       Finally, specify a general equilibrium relaxation and verify that
*       it is calibrated when inialized at the solution of the integrated model:

PARAMETER 
        aref(i,k) Reference output for Armington sector
        dref(i,k,kk) Reference demand
        pdref(i,k) Reference price;

*       The following is an MPSGE model in which we evaluate factor
*       price equilibria given industry structure in each of the
*       individual industries.

$ontext
$model:MCHFGE
$sectors:
        D(i,k) ! Aggregate demand (CES)
        Y(i,k) ! Industry supply

$commodities:
        P(i,k) ! Composite price
        C(i,k) ! Output price
        W(f,k) ! Factor wages

$consumers:
        R(k) ! Representative agent

*       Agents are endowed with factors and demand goods. This
*       demand function is exogenously specified and unchanged
*       through the iterative procedure:

$demand:R(k) s:1
        e:W(f,k)        q:s(f,k)
        d:P(i,k)        q:alpha(i)

*       Producers combine primary factor inputs to produce
*       commodity outputs. These cost functions are exogenous
*       and unchanged in the iterative sequence.

$prod:Y(i,k) S:1
        o:C(i,k)        Q:1
        i:W(f,k)        q:beta(i,f) p:1

*       Iterative adjustments of the GE model occur within
*       the "implicit Armington aggregator".

*       Domestic demand substitutes across goods from
*       different firms. AREF, DREF and PDREF are based
*       on current equilibrium values from the industry-
*       level partial equilibrium models. Increases in
*       the total number of firms supplying market k
*       are reflected in the value of aref(). Changes in
*       the sourcing of inputs from different regions are
*       represented by the reference price-quantity pairs.
*       The underlying model is based on firm-level
*       product differentiation, so while this model appears
*       to have the structure of a conventional Armington
*       model, entry and exit decisions will endogenously
*       alter the implied structure of the Armington model.

$prod:D(i,k) s:sigma
        o:P(i,k)        q:aref(i,k)
        i:C(i,kk)       q:dref(i,kk,k) p:pdref(i,kk)
$offtext
$sysinclude mpsgeset MCHFGE

*       We now verify that we can compute equilibria using
*       block Gauss-Seidel decomposition. These calculations
*       involve partial equilibrium models (one for each commodity)
*       followed by a general equilibrium calculation which assesses
*       the comparative advantage effects in a seemingly conventional
*       CRTS Armington trade model.
*       In this program we first compute the equilibrium with
*       this decomposition procedure, and we then evaluate precision
*       of the estimation in an integrated model.

PT.L(i,k,kk) = 1;
QT.L(i,k,kk) = 1;
N.L(i,k,kk) = 1;
PHIS.L(i,k,kk) = 1;
PHIT.L(i,k,kk) = 1;
R.L(k) = sum(f,s(f,k));
Y.L(i,k) = alpha(i) * R.L(k);
M.L(i,k) = Y.L(i,k);
C.LO(i,k) = 0.00001;
W.LO(f,k) = 0.00001;
P.LO(i,k) = 0.00001;
PT.LO(i,k,kk) = 0.00001;
PHIS.LO(i,k,kk) = b;
PHIT.LO(i,k,kk) = b;

*       Initialize productivity levels using autarky values:

PHIS.L(i,k,kk) = 0;
PHIT.L(i,k,kk) = 0;
PHIS.L(i,k,kk)$t(i,k,kk) = b * (fc(i,k,kk)/(fe(i)*delta) * (sigma-1)/(a+1-sigma))**(1/a);
PHIT.L(i,k,kk)$t(i,k,kk) = (a/(a+1-sigma))**(1/(sigma-1)) * PHIS.l(i,k,kk);
THETA.L(i,k,kk)$t(i,k,kk) = 1;

*       The partial equilibrium relaxation has fixed income
*       levels and fixed numbers of firms:

alias (i,ii);

*       Partial equilibrium solution holds marginal cost of commodities, income
*       and aggregate production of each commodity fixed:

R.FX(k) = R.L(k);
Y.L(i,k) = alpha(i) * R.L(k)/C.L(i,k);

set iter Outer iterations /iter0*iter5/;

parameter       iterlog         Iteration log
                nfirm           Iteration log (number of firms);

loop(iter,

*       Produce an iteration log:

        iterlog(iter,"R"," ",k) = R.L(k);
        iterlog(iter,"Y",i,k) = Y.L(i,k);
        iterlog(iter,"M",i,k) = M.L(i,k);
        iterlog(iter,"C",i,k) = C.L(i,k);

*       One subproblem for each commodity:
        loop(ii,
          eq(i) = no;
          eq(ii) = yes;
          t(i,k,kk) = no;
          t(ii,k,kk) = yes$(1/tau(ii,k,kk));
          SOLVE MCHFPE using MCP;
        );

        nfirm(iter,i) = M.L(i,"h");

*       Recalibrate the general equilibrium model to the
*       partial equilibrium solution:

        aref(i,k) = alpha(i)*R.L(k)/P.L(i,k);
        dref(i,k,kk) = PT.L(i,k,kk)*QT.L(i,k,kk)*N.L(i,k,kk)/C.L(i,k);
        pdref(i,k) = C.L(i,k);

*       Remove bounds to compute the general equilibrium solution:

        Y.l(i,k) = sum(kk, dref(i,k,kk));
        R.UP(k) = +inf; R.LO(k) = 0;
        R.FX(k)$(ord(k)=1) = R.L(k);

$include MCHFGE.GEN
        SOLVE MCHFGE USING MCP;

*       Fix prices and activity levels which enter the partial equilibrium
*       model as exogenous parameters:

        t(i,k,kk) = yes$(1/tau(i,k,kk));
        R.FX(k) = R.L(k);
);

*       Check this solution using the fully-specified model:

eq(i) = yes;
t(i,k,kk) = yes$(1/tau(i,k,kk));
R.UP(k) = +inf; R.LO(k) = 0;
MCHF.ITERLIM = 0;
SOLVE MCHF USING MCP;

*       Generate a plot:

$libinclude plot nfirm

*       Display iterative sequence:

option iterlog:1:1:3;
display iterlog;


----    802 PARAMETER iterlog  Iteration log

             R..H        R..F       Y.1.H       Y.1.F       Y.2.H       Y.2.F       M.1.H       M.1.F       M.2.H

iter0      2200.0      2200.0      1100.0      1100.0      1100.0      1100.0      1100.0      1100.0      1100.0
iter1      2200.0      2200.0      1130.4      1061.4      1061.4      1130.4      4767.8      4767.8      4767.8
iter2      2200.0      2200.0      1146.6      1045.7      1045.7      1146.6      4955.1      4546.7      4546.7
iter3      2200.0      2200.0      1148.2      1044.2      1044.2      1148.2      4975.2      4527.3      4527.3
iter4      2200.0      2200.0      1148.3      1044.1      1044.1      1148.3      4977.1      4525.5      4525.5
iter5      2200.0      2200.0      1148.3      1044.0      1044.0      1148.3      4977.3      4525.3      4525.3

    +       M.2.F       C.1.H       C.1.F       C.2.H       C.2.F

iter0      1100.0         1.0         1.0         1.0         1.0
iter1      4767.8         1.0         1.0         1.0         1.0
iter2      4955.1         1.0         1.0         1.0         1.0
iter3      4975.2         1.0         1.0         1.0         1.0
iter4      4977.1         1.0         1.0         1.0         1.0
iter5      4977.3         1.0         1.0         1.0         1.0
>