Firm dynamics with an example of exporter dynamics

1 What is firm dynamics?

The firm dynamics literature wants to study firms’ behavior over time. Dynamic models were developed to study firms’ entry and exit decisions, age, size distribution over time and even R& D, export and FDI decisions, etc. More recently, spatial components are incorporated to study, say, firms’ heterogeneous location choices for their subsidiaries.

To dig into the literature, a good starting point is Hopenhayn (1992), which provides a general theoretical framework that features uncertainty towards future profits, costly market entry and endogenous exit. The literature then expands dramatically with various focuses. Just to list a few, Klette and Kortum on innovations and firm distribution, Cooley and Quadrini (2001) on financial frictions, Das et al. (2007) on exporter dynamics.

2. What’s specific about exporter dynamics?

In a world where large multinationals take up a huge chunk of global economy via exports, the roles these giants play is a very interesting research topic. The OECD estimated in 2018 that multinational corporations account for half of global exports, nearly a third of world GDP (28%), and about a fourth of global employment. They hold great market power and sometimes even directly alter domestic and foreign policies. Understanding their growth and policy implications is a fascinating ongoing topic, which requires a lot of the tools discussed here. We can use the dynamic trade framework to study firms’ individual and aggregate behavior, international policies and their implications on consumer welfare.

For the remaining part of this section, I will discuss some empirical evidence about firm-level exporting behavior. I mainly follow from the literature review on firm dynamics and trade by George, Costas and Kim, ARE 2021.

3. Firm-level facts (from firm-level data of US and Columbia):

Static or time-invariant facts: * There are few exporters. * Exports are a small share of an exporter’s total sales. * Exporters are relatively large firms.

These facts can already trigger our thinking about the underlying firm distribution in size and productivity. These are consistent with the results of Melitz (2003). Within a group of firms originally producing domestically, only the most productive firms enter the export market. They stay in the market and force the less productive firms to exit.

If we further take a look at the dynamic aspect of the data, there are some more interesting observations:

• Past export participation is the main predictor of current export participation.
• Exporter are less likely to exit if they are older and larger exporters in the past.
• Entry rate is increasing in size and past export activity.
• The longer an exporter stays in the market, the more it tends to export.

These can explain why the dynamic aspect is vital to studying the exporting firms’ behavior. I will show some predictions of the results when we go to the theoretical framework.

There are also evidences about heterogeneity between exporters on their export destinations, shipment volumes and frequencies and inventories.

At last, The long-run response of aggregate trade volumes to changes in trade policy is larger than the short-run response. An important reason, as shown in Kehoe and Ruhl (2013), is that the aggregate trade growth is stronger in the long term due to the extensive margin. That is, during a trade liberalization, previously non-traded products grow faster than products that have already been traded between the same country pairs. One take-away is that larger quantities of firms (for concreteness, plants) would become exporters when there are ease in bilateral trade policies. Combining with the fact that exporters take time to grow, this gives a possible explanation for the long-run growth of trade.

4. Some questions to think about:

On building a model, what structures do we need on individual firm’s problem to retrieve the same results as we observe? For example, how to make sure that only a few firms enter into the export market and coincidentally they are relatively large? What about on the market structure and the demand side?

Next, we can think about idiosyncratic and aggregate uncertainties. Idiosyncratic shock is convenient to generate endogenous entry and exit, as in the firm dynamics literature. Now we can look at if there are any heterogeneous short and long run responses to aggregate uncertainties, say real exchange rate fluctuations. If so, what are the policy implications?

A very important topic in the trade literature is the gains from trade. That is the welfare gain of household in a country facing some sort of trade liberalization. Although the focus of our model is mainly on the firm side, we can still analyze the aggregate impact under a general equilibrium framework.

5. A canonical model for exporter dynamics

The model is a slightly tweaked version of Das et al. (2007), based upon Melitz (2003). For further reference, Ruhl and Willis (2017) and Alessandria et al. (2021).

This section is arranged as follows. I will first lay out the static model and provide analytical solutions to it. Next, I will provide the dynamic setup of the model and compute the model. To look at the firms’ short-term reaction to some sudden and unexpected (MIT) shocks, I will then illustrate the algorithms for obtaining the transition dynamics. In the end, I will provide some basic calibration and some preliminary results. I run some experiments and look at the exporters’ full dynamics facing different shocks.

5.1 Static model

5.1.1 Consumer

A representative consumer in domestic economy has preferences over an aggregate consumption good, \(C\). Each index \(j\) denotes a differentiated variety. \(C\) takes the form: \[ C = \left(\sum_{j=1}^{J} c_j^\frac{\theta - 1}{\theta} \right)^\frac{\theta}{\theta - 1}, \] where \(\theta\) is the elasticity of substitution between varieties, and J is the number of available varieties. The consumer maximize his utility subject to the budget constraint: \[ \sum_{j=1}^{J} c_j p_j = I. \] The consumer takes prices as given and chooses optimal consumption \(c_j\) of each variety: \[ c_j = (\frac{p_j}{P})^{-\theta}C, \] where P is an aggregate price level given by: \[ P = \left(\sum_{j=1}^{J} p_j^{1-\theta} \right)^{\frac{1}{1 - \theta}}. \] The foreign market is identical to the domestic market with demand and consumption satisfying: \[ \begin{aligned} c_j^* &= (\frac{p_j^*}{P^*})^{-\theta}C^*, \\ P^* &= \left(\sum_{j=1}^{J} {p_j^*}^{1-\theta} \right)^{\frac{1}{1 - \theta}}. \\ \end{aligned} \] Countries have same elasticity \(\theta\) for differentiated goods but differ in aggregate level of demands and prices.

5.1.2 Firm’s static problem

Assume each firm operates solely on a market \(j\), so “plants” will be a better term than “firms”. Plant operates monopolistically in each differentiated variety. Plant chooses the optimal amounts to produce in domestic and foreign market. A unit good produced with plant-level technology $_j $, labor \(n_j\) and capital \(k_j\) subject to: \[ f(\tilde{\epsilon}_j, n_j, k_j) = \tilde{\epsilon}_j n_j^{\alpha_N}k_j^{\alpha_K}, \] where we assume nonincreasing returns to scale of the production technology, \(\alpha_N + \alpha_K <= 1\).

In each period, the plant chooses the prices, production, inputs, and export status next period \(X'\) (\(X'=1\) if exporting next period and \(X'=0\) if not) to maximize its infinite horizon utility. The plant also faces an iceberg cost \(\xi_j\), where fraction \(\xi_j - 1\) of an export shipment is destroyed in transportation. In addition, foreign market might put an ad valorem tariff \(\tau_j\) on exports.

We can thus first solve for the plant’s single period or static problem given all the states and plant’s export decision last period. A plant’s profit is sum of revenue from domestic and foreign sales less the input costs: \[ \Pi_j = p_j y_j + I(X_j = 1) (1 - \tau_j) Q p_j^* y_j^* - wn_j - rk_j, \] where \(Q\) is the real exchange rate; r is the rental rate of capital; and w is the wage. The plant is also subject to its feasibility constraint: \[ y_j + \xi_j y_j^* = \tilde{\epsilon}_j n_j^{\alpha_N}k_j^{\alpha_K}. \] With market clearing condition in both markets \(c_j = y_j\) and \(c^*_j = y^*_j\) if \(X_j = 1\), with demand functions from the household problem. The plant charges constant markup over marginal cost in both markets: \[ \begin{aligned} p_j &= \frac{\theta}{\theta - 1} {MC}_j \\ Q p_j^* &= \frac{1}{1- \tau_j} \frac{\theta}{\theta - 1} \xi_j {MC}_j. \end{aligned} \] where \(MC_j = \frac{1}{\tilde{\epsilon_j}} \left(\frac{w}{\alpha_n}\right)^{\alpha_n} \left(\frac{r}{1 - \alpha_n}\right)^{1 - \alpha_n}\) with \(\alpha_n + \alpha_k = 1\) (if not, the formula for marginal cost will be more complicated).

The plant’s static profit maximization problem is: \[ \Pi_j = \max_{y_j, y_j^*}\{ P C^{\frac{1}{\theta}} y_j^\frac{\theta - 1}{\theta} + I(X_j = 1) (1 - \tau_j) Q P^* {C^*}^{\frac{1}{\theta}} {y_j^*}^\frac{\theta - 1}{\theta} - wn_j - rk_j\}, \] subject to equation the feasibility constraint. The optimal amounts for a plant to produce in domestic and foreign markets are: \[ \begin{aligned} y_{j} &=\frac{ \xi_j^{\theta} (1 - \tau_j)^{-\theta} (\frac{P}{Q P^*})^\theta \frac{C}{C^{*}}}{I \left(X_{j}=1\right) \xi_j + \xi_j^{\theta} (1 - \tau_j)^{-\theta} (\frac{P}{Q P^*})^\theta \frac{C}{C^{*}}} \tilde{\epsilon}_{j} n_{j}^{\alpha_{N}} k_{j}^{\alpha_{K}} \\ y_{j}^{*} &= \frac{I(X_j = 1)}{\xi_j+ \xi_j^{\theta} (1 - \tau_j)^{-\theta} (\frac{P}{Q P^*})^\theta \frac{C}{C^{*}}} \tilde{\epsilon}_{j} n_{j}^{\alpha_{N}} k_{j}^{\alpha_{K}} . \end{aligned} \]

Substituting in \(y_j\) and \(y_j^*\) to the profit maximization problem yields: \[ \begin{aligned} \Pi\left(X_{j}, \xi_j, \epsilon_{j}, Q\right)=&\max _{n_{j}, k_{j}} M n_{j}^{\alpha_{N} \frac{(\theta-1)}{\theta}} k_{j}^{\alpha_{K} \frac{(\theta-1)}{\theta}} - wn_{j}- r k_{j}, \end{aligned} \] where \(M = \left(1+I\left(X_{j}=1\right) \xi_j^{1 - \theta} (1 - \tau_j)^\theta \left(\frac{Q P^*}{P}\right)^{\theta} \frac{C^{*}}{C}\right)^{\frac{1}{\theta}} P C^{\frac{1}{\theta}} \tilde{\epsilon}_{j}^{\frac{\theta-1}{\theta}}\)

From first order conditions of labor and capital, for simplicity write \(a = \alpha_{N} \frac{(\theta-1)}{\theta}\) and \(b = \alpha_{K} \frac{(\theta-1)}{\theta}\), the maximized profit can be solved as: \[ \Pi\left(X_{j}, \xi_j, \epsilon_{j}, Q\right) = M^{\frac{1}{1-a-b}} (1-a-b) (\frac{a}{w})^{\frac{a}{1-a-b}} (\frac{b}{r})^{\frac{b}{1-a-b}}. \] Equate \(\xi_j = 1\) you get exactly the plant’s profit without export iceberg cost. Then with only export fixed cost driving the export barrier, the model is equivalent to the baseline sunk cost model from . Assuming no tariff, using the same notation, the profit a plant earns with an iceberg cost is exactly \(\left({\frac{1+I\left(X_{j}=1 \right) \xi_j^{1 - \theta} \left(\frac{Q P^*}{P}\right)^{\theta} \frac{C^{*}}{C}}{1+I\left(X_{j}=1\right) \left(\frac{Q P^*}{P}\right)^{\theta} \frac{C^{*}}{C}}}\right)^\frac{1}{\theta (1- a-b)}\) of the same plant without an iceberg cost.

Just for clarification, we can also solve for domestic and foreign sales: \[ \begin{aligned} p_j y_{j} &= \left(\frac{ \xi_j^{\theta} (1 - \tau_j)^{-\theta} (\frac{P}{Q P^*})^\theta \frac{C}{C^{*}}}{I \left(X_{j}=1\right) \xi_j + \xi_j^{\theta} (1 - \tau_j)^{-\theta} (\frac{P}{Q P^*})^\theta \frac{C}{C^{*}}}\right)^{\frac{\theta - 1}{\theta}} \\ & \ \ \ \ \ \ \ \ \ \ P C^{\frac{1}{\theta}} \epsilon_j M^{\frac{a + b}{1-a-b}} (\frac{a}{w})^{\frac{a}{1-a-b}} (\frac{b}{r})^{\frac{b}{1-a-b}} \\ Q p_j^* y_j^* &= I(X_j = 1) \left(\frac{1}{\xi_j+ \xi_j^{\theta} (1 - \tau_j)^{-\theta} (\frac{P}{Q P^*})^\theta \frac{C}{C^{*}}}\right)^{\frac{\theta - 1}{\theta}} \\ & \ \ \ \ \ \ \ \ \ \ Q P^* {C^*}^{\frac{1}{\theta}} \epsilon_j M^{\frac{a + b}{1-a-b}} (\frac{a}{w})^{\frac{a}{1-a-b}} (\frac{b}{r})^{\frac{b}{1-a-b}} . \end{aligned} \]

We normalize the size of domestic aggregate demand \(C\) to 1, so \(C^*\) denotes the size of foreign demand relative the domestic demand. Define \(\epsilon_j = {\tilde{\epsilon}_j}^\frac{\theta - 1}{\theta}\), which follows an idiosyncratic AR1 shock for each plant. We normalize the mean of the \(\epsilon\) process to one.

Firms face idiosyncratic productivity shocks, a fixed cost for becoming an exporter abd and an iceberg cost for exports.

5.2 Firm’s Dynamic problem

We first introduce an export fixed cost structure as in Ruhl and Willis (2017). The plant faces a sunk entry cost \(f_E\) and a continuation cost \(f_C\), \(f_C \leq f_E\). All export decisions are made at the end of previous period and fixed costs are paid upon decision. The exporting cost of each period can be expressed as: \[ f(X, X') = (1-X) * f_E * X' + X * f_C * X' . \]

We then link the export decision of a plant with an endogenous and stochastic iceberg transportation cost \(\xi_j\) as in Alessandria et al. (2021). A non-exporting plant doesn’t ship any goods abroad and we set its iceberg cost this period to infinity. The nonexporter can deterministically lower its fixed cost next period to \(\xi_E\) by paying \(f_E\). If a current exporter pays \(f_C\) to continue exporting next period, it draws iceberg cost \(\xi_j \in \{\xi_E, \xi_C\}\) next period stochastically. The structure of the iceberg cost conditional on continuing to exporting can be summarized by a markov chain with states \(\xi_E\) and \(\xi_C\) with transition probability \([ \rho_{EE}, 1- \rho_{EE} ; 1 - \rho_{CC}, \rho_{CC} ]\). Assume \(\xi_C \leq \xi_E < \infty\). If an exporter stops paying \(f_C\), it exits the export market and becomes a nonexporter with infinite iceberg cost next period.

We also have a persistent condition as \(\rho_\xi(\xi_E|\xi_C) \leq \rho_\xi(\xi_C|\xi_E)\). With a lot of plants, this particular structure of iceberg cost creates incentives for some plants to invest in export technology and make the “better” exporters to become less susceptible to unfavorable shocks.

One last component is the technological shock \(\epsilon_j\) that governs plant productivity, which generates heterogeneous plants and varying entry and exit export decisions of plants. It follows an idiosyncratic and time-invariant AR(1) processes for each plant: \[\begin{equation} \ln{\epsilon_t} = \rho_\epsilon\ln{\epsilon_{t-1}} + w_{\epsilon,t}, \hspace{1cm} w_\epsilon \sim N(0, \sigma^2_\epsilon). \end{equation}\] The shocks are discretized using Tauchen’s method. There are three key parameters to determine here: {{\(N_\epsilon\), \(\rho_\epsilon\), \(\sigma_\epsilon\)}}. \(N_\epsilon\) is the states of plant’s productivity, currently set to 1000 for simplicity.

We put aside our discussion of exchange rate shocks for the moment as an aggregate shock will make solving for the equilibrium prices and transition dynamics much harder. Exchange rates directly shift the export prices and create more uncertainties which make the foreign market more volatile. In particularly “bad” periods, there can be no firms deciding to export, leading to potentially poor simulation results.

In each period, every plant makes an discrete choice based on their current status and discounted expectation of future profits. The discount rate \(R\) is set to be \(\frac{1}{1+r}\). We characterize the dynamic problem of each plant by the following Bellman equation: \[ V\left(X_{j}, \xi_j, \epsilon_{j}, Q\right)=\max _{X_{j}^{\prime}}\left\{\Pi\left(X_{j}, \xi_j, \epsilon_{j}, Q\right)-f\left(X_{j}, X'_{j} \right)+R \underset{\xi_j^{\prime}, \epsilon_{j}^{\prime}, Q^{\prime}}{\mathbb{E}} V\left(X_{j}^{\prime}, \xi_j^{\prime}, \epsilon_{j}^{\prime}, Q^{\prime}\right)\right\}. \]

The policy function for export entry is thus: \[ X_{j}^{\prime}\left(0, \infty, \epsilon_{j}, Q\right)= \begin{cases}1 & \text{ if } \hspace{0.5cm} - f_E + \\ {} & R \underset{\xi_j^{\prime}, \epsilon_{j}^{\prime}, Q^{\prime}}{\mathbb{E}} \left[V\left(1, \xi_E, \epsilon_{j}^{\prime}, Q^{\prime}\right)-V\left(0, \infty, \epsilon_{j}^{\prime}, Q^{\prime}\right)\right] \geq 0 \\ 0 & \text { otherwise. }\end{cases} \]

The inequality on the right is the condition of a non-exporting plant choosing to be an exporter next period. The first two terms denote the current period profit minus the entry fixed cost. The discounted expectation is the difference of expected future values between starting to become an exporter with a low export technology \(\xi_E\) and staying as an nonexporter.

6. Solving the model

As always, let’s input some parameter value. I made the simplifying assumption that the real exchange rate \(Q\) is fixed.

using LinearAlgebra, Statistics, Distributions
using Plots

Base.@kwdef struct Exporter{T1, T2, R1, S}
    r::R1 = 0.109
    αn::R1 = 0.45

    αk::R1 = 1.0 - αn
    θ::R1 = 5.0
    w::R1 = 0.02

    ξE::R1 = 1.6
    ξC::R1 = 1.2
    ξ::T1 = [ξE; ξC]
    
    ρEE::R1 = 0.92
    ρCC::R1 = 0.92
    Ξ::T2 = [ρEE (1.0-ρEE); (1.0-ρCC) ρCC]

    Q::R1 = 1.0 # fix exchange rate for the moment

    fE::R1 = 0.7
    fC::R1 = 0.35
    
    ρϵ::R1 = 0.872524
    σϵ::R1 = 0.115886

    Cstar::R1 = 0.7 #C normalized to be 1.0

    grids::S = 100

    R = (1/(1+r))^0.25 # quarterly discount factor
end

E = Exporter() 

Exporter{Vector{Float64}, Matrix{Float64}, Float64, Int64}(0.109, 0.45, 0.55, 5.0, 0.02, 1.6, 1.2, [1.6, 1.2], 0.92, 0.92, [0.92 0.07999999999999996; 0.07999999999999996 0.92], 1.0, 0.7, 0.35, 0.872524, 0.115886, 0.7, 100, 0.9744669483879411)

6.1 Discretize an AR(1) process

Recall that the reason we use discrete approximations of value functions is that the computer can’t process continuos variables. Same logic applies with an AR(1) process. We want to approximate the AR(1) process with a discrete markov chain, which is characterized by discrete states and a transition probability matrix. Check Quantecon’s lecture notes for a detailed explanation.

Again the AR(1) process we want to discretize is: \[ \ln{\epsilon_t} = \rho_\epsilon\ln{\epsilon_{t-1}} + w_{\epsilon,t}, \hspace{1cm} w_\epsilon \sim N(0, \sigma^2_\epsilon). \]

The code (adopted from Quantecon) is:

function tauchen(N::Integer, ρ::T1, σ::T2, μ=zero(promote_type(T1, T2)), n_std::T3=3) where {T1 <: Real, T2 <: Real, T3 <: Real}
  # Get discretized space
  a_bar = n_std * sqrt(σ^2 / (1 - ρ^2))
  y = range(-a_bar, stop=a_bar, length=N)
  d = y[2] - y[1]

  # Get transition probabilities
  Π = zeros(promote_type(T1, T2), N, N)
  for row = 1:N
      # Do end points first
      Π[row, 1] = cdf(Normal(),(y[1] - ρ*y[row] + d/2) / σ)
      Π[row, N] = 1.0 - cdf(Normal(),(y[N] - ρ*y[row] - d/2) / σ)

      # fill in the middle columns
      for col = 2:N-1
          Π[row, col] = (cdf(Normal(),(y[col] - ρ*y[row] + d/2) / σ) -
                        cdf(Normal(),(y[col] - ρ*y[row] - d/2) / σ))
      end
  end

  yy = exp.(y .+ μ / (1 - ρ)) # center process around its mean (wbar / (1 - rho)) in new variable
  # take exponential to get to a ln process

  # renormalize. In some test cases the rows sum to something that is 2e-15
  # away from 1.0, which caused problems in the MarkovChain constructor
  Π = Π./sum(Π, dims = 2)

  return Π, yy
end

tauchen (generic function with 3 methods)

Π, ϵ_vec = tauchen(100, E.ρϵ, E.σϵ)

([0.2355397082043846 0.0398220365205688 … 0.0 0.0; 0.2035639923092491 0.03686888293414207 … 0.0 0.0; … ; 4.784202347164862e-30 1.4778853656337744e-29 … 0.03686888293414202 0.20356399230924901; 1.3827115757372961e-30 4.347177100735715e-30 … 0.039822036520568704 0.23553970820438463], [0.4908675431085489, 0.4979749102029403, 0.5051851862545925, 0.5124998613024105, 0.5199204469598645, 0.5274484767273707, 0.5350855063091983, 0.5428331139349621, 0.5506929006857724, 0.5586664908251056  …  1.7899766970506505, 1.8158941194896647, 1.84218680535361, 1.8688601881549594, 1.8959197800790755, 1.9233711731233287, 1.9512200402527142, 1.9794721365721941, 2.0081333005160213, 2.0372094550542794])

6.2 Functions

Exporter fixed cost (just for illustration, we won’t use this function)

function expfixcost(X1::Int, X2::Int; E = E)
    #X1 and X2 are 0, 1 variables
    f = (1-X1) * E.fE * X2 + X1 * E.fC * X2;
    return f
end

expfixcost (generic function with 1 method)

Firm’s profit function

function profit(X, Evalue, ξ_idx, P, Pstar; E= E, Qvalue = E.Q, τ = 0.0)
    (;θ, Cstar, αn, αk, w, r, ξ ) = E
    M = (1.0 + X * ξ[ξ_idx]^(1.0 - θ)* (1.0 - τ)^θ * ((Qvalue * Pstar/ P)^θ) * Cstar)^(1.0/θ) * P * Evalue
    a = αn * (θ - 1.0)/θ
    b = αk * (θ - 1.0)/θ
    Profit = (1.0-a-b) * M^(1.0/(1.0-a-b)) * (a/w)^(a/(1.0-a-b)) * (b/r)^(b/(1.0-a-b))
    
    return Profit
end

profit (generic function with 1 method)

Sales at domestic and foreign markets:

function sales(X, Evalue, ξ_idx, P, Pstar; E= E, Qvalue = E.Q, τ = 0.0)
    (;θ, Cstar, αn, αk, w, r, ξ) = E
    M = (1.0 + X * ξ[ξ_idx]^(1.0 - θ)* (1.0 - τ)^θ * ((Qvalue * Pstar/ P)^θ) * Cstar)^(1.0/θ) * P * Evalue
    a = αn * (θ - 1.0)/θ
    b = αk * (θ - 1.0)/θ
    con = ξ[ξ_idx]^θ * (1.0 - τ)^(-θ) * (P/ (Qvalue * Pstar))^θ/Cstar
    sales_do = (con/(X * ξ[ξ_idx] + con))^((θ-1.0)/θ) * P * Evalue * M^((a + b)/(1.0-a-b)) * (a/w)^(a/(1.0-a-b)) * (b/r)^(b/(1.0-a-b))
    if X == 0
        sales_fo = 0.0
    else
        sales_fo = (1.0-τ) * (1.0/(ξ[ξ_idx] + con))^((θ-1.0)/θ) * Qvalue * Pstar * Cstar^(1.0/θ) * Evalue * M^((a + b)/(1.0-a-b)) * (a/w)^(a/(1.0-a-b)) * (b/r)^(b/(1.0-a-b)) # foreign sales net of tariffs
    end
    return [sales_do, sales_fo]
end

sales (generic function with 1 method)

Take a try on the functions. Again . is used for broadcasting and Ref() is used to fix inputs.

profit.([1, 0], [1.2, 0.8], [1, 1], Ref(0.1), Ref(0.2))

2-element Vector{Float64}:
 0.0860809333951814
 0.0025658309095826434

reduce(hcat, sales.([1, 0], [1.2, 0.8], [1, 1], Ref(1.0), Ref(1)))'

2×2 adjoint(::Matrix{Float64}) with eltype Float64:
 9742.14  1040.57
 1282.92     0.0

Individual prices (\(p_j\) and \(p_j^*\))

The function defined here needs to be modified to incorporate cases where \(\alpha_n + \alpha_k \neq 1\).

function ind_price(X, Evalue, ξ_idx; E = E, Qvalue = E.Q, τ = 0.0)
    (;θ, αn, αk, w, r, ξ) = E
    MC = Evalue^(θ/(1.0 - θ)) * (w/αn)^αn * (r/αk)^αk
    do_p = (θ/(θ - 1.0)) * MC
    if X == 0
        fo_p = -1.0 # make foreign price negative for non-exporters
    else
        fo_p = (θ/(θ - 1.0)) / (1.0 - τ) * ξ[ξ_idx] * MC/Qvalue
    end
    return [do_p, fo_p]
end

ind_price (generic function with 1 method)

ind_p_mat = reduce(hcat, ind_price.([0 1], [1.0 1.2], [1 1]))'[:,1]
ind_pstar_mat = reduce(hcat, ind_price.([0 1], [1.0 1.2], [1 1]))'[:,2]

2-element Vector{Float64}:
 -1.0
  0.16104839447304678

Aggregate price \(P\) and \(P^*\)

function agg_price(ind_p_mat, ind_pstar_mat; θ = E.θ)
    P = (sum(ind_p_mat.^(1.0 - θ)))^(1.0/(1.0 - θ))
    Pstar = ((sum(ind_pstar_mat[ind_pstar_mat .>= 0]).^(1.0 - θ)))^(1.0/(1.0 - θ))
    return P, Pstar
end

agg_price (generic function with 1 method)

agg_price(ind_p_mat, ind_pstar_mat)

(0.09250365544979017, 0.16104839447304678)

6.3 Solve the firm’s problem

Solving the model is rather standard except that I chose to separate out the problem of a non-exporter and an exporter with either exporting technology \(\xi\).

function solve_exporter_fixQ(E, P, Pstar, τ; vNX0 = zeros(E.grids, 1), vEX0 = zeros(E.grids, length(E.ξ)), tol = 1e-6, Q = E.Q)
    # unpacking
    (; R, fE, fC, ξ, Ξ, fE, fC, ρϵ, σϵ, grids) = E

    vNX1 = similar(vNX0) # to store value of non-exporters
    vEX1 = similar(vEX0) # to store value of exporters
    polNX = similar(vNX0, Int) # to store policy of non-exporters
    polEX = similar(vEX0, Int) # to store policy of exporters

    REvNX = similar(vNX0) # expected value if choosing not to export
    REvEX_NX = similar(vNX0) # expected value if choosing to export as a non-exporter
    REvEX_EX = similar(vEX0) # expected value if choosing to export as an exporter with different ξ's
    
    domesticm, ~ = sales(0, 1.0, 1, P, Pstar; E = E, Qvalue = Q, τ = 0.0) # sales of a non-exporting medium firm

    # productivity process
    Π, ϵ_vec = tauchen(grids, ρϵ, σϵ)

    # Precompute the profits for firms at all states 
    proNX = profit.(zeros(Int, grids, 1), ϵ_vec, ones(Int, grids, 1), Ref(P), Ref(Pstar); E= E, Qvalue = Q, τ = τ)
    # input for ξ doesn't matter

    proEX = similar(vEX0)
    for i in 1:length(ξ)
        proEX[:,i] = profit.(ones(Int, grids, 1), ϵ_vec, i.*ones(Int, grids, 1), Ref(P), Ref(Pstar); E= E, Qvalue = Q, τ = τ)
    end

    # VFI
    iter = 0
    while true
        distance = zero(eltype(vNX0))
        iter += 1 

        # Precompute RE[V(.)|X, ϵ, ξ] for better performance
        REvNX = R * Π * vNX0
        REvEX_NX = R * Π * vEX0[:,1] # enters with entry export technology
        REvEX_EX = R * Π * vEX0 * Ξ'

        vNX1 = max.((proNX + REvNX), (proNX .- fE * domesticm + REvEX_NX))
        vEX1 = max.((proEX .+ REvNX), (proEX .- fC * domesticm + REvEX_EX))

        distance = max(maximum(abs.(vNX1 - vNX0)), maximum(abs.(vEX1 - vEX0)))

        (distance < tol || iter == 2000) && break
        vNX0 .= vNX1
        vEX0 .= vEX1
    end

    # get policy
    polNX .= 0
    polEX .= 0
    for i in 1:grids
        if (REvNX[i]) < (- fE * domesticm + REvEX_NX[i])
            polNX[i] = 1
        end
        for j in 1:length(ξ)
            if (REvNX[i]) < (- fC * domesticm + REvEX_EX[i,j])
                polEX[i,j] = 1
            end
        end
    end

    return (vNX = vNX1, vEX = vEX1, polNX = polNX, polEX = polEX, Π = Π, ϵ_vec =ϵ_vec, iter = iter)

end

solve_exporter_fixQ (generic function with 1 method)

@time begin
    sol_try = solve_exporter_fixQ(E, 0.1, 0.1, 0.0);
end

  0.230753 seconds (373.87 k allocations: 31.695 MiB, 3.40% gc time, 65.36% compilation time)

(vNX = [0.43476592047069756; 0.43593320133695335; … ; 1.4557350993346505; 1.4970252086816265;;], vEX = [0.43478975561531463 0.4348412512981293; 0.43595881274862097 0.4360141460454338; … ; 1.5042340280033728 1.641819695915486; 1.547864034400909 1.6930444062885848], polNX = [0; 0; … ; 1; 1;;], polEX = [0 0; 0 0; … ; 1 1; 1 1], Π = [0.2355397082043846 0.0398220365205688 … 0.0 0.0; 0.2035639923092491 0.03686888293414207 … 0.0 0.0; … ; 4.784202347164862e-30 1.4778853656337744e-29 … 0.03686888293414202 0.20356399230924901; 1.3827115757372961e-30 4.347177100735715e-30 … 0.039822036520568704 0.23553970820438463], ϵ_vec = [0.4908675431085489, 0.4979749102029403, 0.5051851862545925, 0.5124998613024105, 0.5199204469598645, 0.5274484767273707, 0.5350855063091983, 0.5428331139349621, 0.5506929006857724, 0.5586664908251056  …  1.7899766970506505, 1.8158941194896647, 1.84218680535361, 1.8688601881549594, 1.8959197800790755, 1.9233711731233287, 1.9512200402527142, 1.9794721365721941, 2.0081333005160213, 2.0372094550542794], iter = 375)

Find the cutoff productivity levels

function find_cutoff_prod(sol)
    (; polNX, polEX, ϵ_vec) = sol
    cutoff_ϵ = zeros(eltype(ϵ_vec), size(polEX)[2] + 1)
    for i in 1:size(polNX)[1]
        if polNX[i] == 1
            cutoff_ϵ[1] = ϵ_vec[i]
            break
        end
    end
    
    for j in 1:size(polEX)[2]
        for i in 1:size(polEX)[1]
            if polEX[i, j] == 1
                cutoff_ϵ[j + 1] = ϵ_vec[i]
                break # only break out of the inner loop
            end
        end
    end

    return cutoff_ϵ
end

find_cutoff_prod (generic function with 1 method)

Now we can find out the cutoff levels for the productivity \(\epsilon\) for different export status \(X\) and export technology \(\xi\) .

find_cutoff_prod(sol_try)

3-element Vector{Float64}:
 2.0081333005160213
 1.6658337302526012
 1.342716960827644

This says that among different states, non-exporter needs the highest draw of productivity to become an exporter and the exporter with highest \(\xi\) will not go back to a non-exporter only if the productivity falls very low.

6.4 Stationary distribution for firms given prices

function compute_stationary_distribution(sol;E = E, tol = 1e-8, pdf_0 = fill(1.0 / (size(sol_try.vEX)[1] * (size(sol_try.vEX)[2] + 1))  , (size(sol_try.vEX)[1] , (size(sol_try.vEX)[2] + 1))))

    (; polEX, polNX, vEX, vNX) = sol

    (; Ξ, ρϵ, σϵ, grids) = E

    # productivity process
    Π, ϵ_vec = tauchen(grids, ρϵ, σϵ)
    
    pdf_1 = similar(pdf_0)

    iter = 0
    while true 
        fill!(pdf_1, zero(eltype(pdf_0)))
        iter +=1

        for i in 1:grids
            for iprime in 1:grids
                pdf_1[iprime,1] += (1 - polNX[i]) * pdf_0[i, 1] * Π[i, iprime]
                pdf_1[iprime,2] += polNX[i] * pdf_0[i, 1] * Π[i, iprime]
                for j in 1:size(Ξ)[1]
                    pdf_1[iprime,1] += (1 - polEX[i,j]) * pdf_0[i, j + 1] * Π[i, iprime]
                    for jprime in 1:size(Ξ)[1]
                        pdf_1[iprime,jprime + 1] += polEX[i, j] * Ξ[j, jprime] * Π[i, iprime] * pdf_0[i ,j + 1] 
                    end
                end
            end
        end
        
        distance = zero(eltype(pdf_0))
        for (a, b) in zip(pdf_0, pdf_1)
            distance = max(abs(a - b), distance)
        end
        
        (distance < tol || iter == 2000) && break 
        pdf_0 .= pdf_1
    end
    pdf_1 = pdf_1
    return pdf_1, iter
end

compute_stationary_distribution (generic function with 1 method)

@time begin
    pdf_stationary, _ = compute_stationary_distribution(sol_try)
end

  0.082796 seconds (104.75 k allocations: 6.769 MiB, 93.08% compilation time)

([0.0011570502171430112 1.134093337951658e-22 1.3041790003368101e-21; 0.0002998023372470313 2.0936653656776906e-22 2.4076453556327584e-21; … ; 0.0002164516680350745 7.504147054101457e-5 8.309198670945468e-6; 0.0007172049628995025 0.00040114335735656076 3.870189688694214e-5], 49)

Check if it is still a distribution:

sum(pdf_stationary)

1.0000000000000007

Look at the shares of exporters at each states.

NX_share, EX_low_share, EX_high_share = sum(pdf_stationary, dims = 1)
println("The share of nonexporters is $NX_share, the share of exporters with low technology is $EX_low_share, the share of exporters with high technology is $EX_high_share, given prices in the stationary distribution, ")

The share of nonexporters is 0.9966134925743302, the share of exporters with low technology is 0.002694855795406125, the share of exporters with high technology is 0.000691651630264524, given prices in the stationary distribution, 

6.5 Stationary prices

Now we have all the pieces ready for the stationary prices.

The algorithm goes as follows:

• Fix the number of firms N
• Start from guesses of \(P_0\) and \(P^*_0\)
• For each \(P_t\) and \(P^*_t\) solve the firm’s problem given prices.
• Compute invariant distribution \(\lambda(X, \epsilon, \xi)\)
• Get new individual firm prices according to first order conditions of firms.
• Get the updated aggregate prices \(P'_t\) and \(P'^{*}_t\).
• Guess \(P_{t + 1} = ϵ P'_t + (1 - ϵ) P_t\) and \(P^*_{t + 1} = ϵ P'^*_t + (1 - ϵ) P^*_t\). Go back to the second step and iterate until convergence.

function compute_stationary_equilibrium(E; P0=0.1, Pstar0=0.1, tol=1e-5, ϵ=0.2, τ=0.0, N=2000)
    (; Ξ, ρϵ, σϵ, grids, θ) = E

    Π, ϵ_vec = tauchen(grids, ρϵ, σϵ) # discretized AR1

    ind_p_mat = zeros(grids, size(Ξ)[1] + 1)
    ind_pstar_mat = zeros(grids, size(Ξ)[1] + 1)

    sol = solve_exporter_fixQ(E, P0, Pstar0, τ) #solve the DP under some guesses of aggregate prices
    pdf, ~ = compute_stationary_distribution(sol) #compute the stationary distribution using the policy

    for i in 1:grids
        ind_p_mat[i, 1], ind_pstar_mat[i, 1] = ind_price(0, ϵ_vec[i], 1; τ=τ) .* pdf[i, 1] #compute individual prices for non-exporting firms
        for j in 1:size(Ξ)[1]
            ind_p_mat[i, j+1], ind_pstar_mat[i, j+1] = ind_price(1, ϵ_vec[i], j; τ=τ) .* pdf[i, j+1] #compute individual prices for exporters
        end
    end

    Pprime1 = (sum(ind_p_mat.^(1.0 - θ) .*pdf .* N ))^(1.0/(1.0 - θ)) # compute aggregate domestic price
    temp = ind_pstar_mat.^(1.0 - θ) .*pdf .* N 
    Pstarprime1 = (sum(temp[temp .> 0.0]))^(1.0/(1.0 - θ)) # compute aggregate foreign price of exporters

    P1 = ϵ * Pprime1 + (1.0 - ϵ) * P0 # update rule for domestic price
    Pstar1 = ϵ * Pstarprime1 + (1.0 - ϵ) * Pstar0 # update rule for foreign price

    P0 = copy(P1)
    Pstar0 = copy(Pstar1)

    iter = 0
    while true #put the steps above in iteration and iterate until convergence
        distance = zero(eltype(P0))
        iter += 1

        sol = solve_exporter_fixQ(E, P0, Pstar0, τ)
        pdf, ~ = compute_stationary_distribution(sol)

        for i in 1:grids
            ind_p_mat[i, 1], ind_pstar_mat[i, 1] = ind_price(0, ϵ_vec[i], 1; τ=τ)
            for j in 1:size(Ξ)[1]
                ind_p_mat[i, j+1], ind_pstar_mat[i, j+1] = ind_price(1, ϵ_vec[i], j; τ=τ)
            end
        end

        Pprime1 = (sum(ind_p_mat.^(1.0 - θ) .*pdf .* N ))^(1.0/(1.0 - θ))
        temp = ind_pstar_mat.^(1.0 - θ) .*pdf .* N 
        Pstarprime1 = (sum(temp[temp .> 0.0]))^(1.0/(1.0 - θ))

        P1 = ϵ * Pprime1 + (1.0 - ϵ) * P0
        Pstar1 = ϵ * Pstarprime1 + (1.0 - ϵ) * Pstar0

        distance = max(abs(P1 - P0), abs(Pstar1 - Pstar0))
        (distance < tol || iter == 2000) && break

        P0 = copy(P1)
        Pstar0 = copy(Pstar1)
    end

    return P1, Pstar1, ind_p_mat, ind_pstar_mat, sol, pdf, iter
end

compute_stationary_equilibrium (generic function with 1 method)

@time begin
    P_rce, Pstar_rce, ind_p_mat_rce, ind_pstar_mat_rce, sol_rce, pdf_rce, iter_rce = compute_stationary_equilibrium(E)
end

  0.531602 seconds (407.32 k allocations: 100.759 MiB, 3.33% gc time, 15.25% compilation time)

(0.015952055070459048, 0.022740007481633726, [0.30768588844713596 0.30768588844713596 0.30768588844713596; 0.30220639279011097 0.30220639279011097 0.30220639279011097; … ; 0.052883757118820714 0.052883757118820714 0.052883757118820714; 0.05194196443888268 0.05194196443888268 0.05194196443888268], [-1.0 0.4922974215154175 0.36922306613656314; -1.0 0.4835302284641776 0.3626476713481332; … ; -1.0 0.08461401139011314 0.06346050854258485; -1.0 0.08310714310221229 0.062330357326659216], (vNX = [2.6372572588229866e-5; 2.6534988291849078e-5; … ; 0.00019073299008637822; 0.0001965218456753796;;], vEX = [2.6387095351844466e-5 2.6418471692987123e-5; 2.6550593335722793e-5 2.6584307936684532e-5; … ; 0.00021410337885572679 0.00030581855386856; 0.00022130259243730753 0.00031748413001791764], polNX = [0; 0; … ; 1; 1;;], polEX = [0 0; 0 0; … ; 1 1; 1 1], Π = [0.2355397082043846 0.0398220365205688 … 0.0 0.0; 0.2035639923092491 0.03686888293414207 … 0.0 0.0; … ; 4.784202347164862e-30 1.4778853656337744e-29 … 0.03686888293414202 0.20356399230924901; 1.3827115757372961e-30 4.347177100735715e-30 … 0.039822036520568704 0.23553970820438463], ϵ_vec = [0.4908675431085489, 0.4979749102029403, 0.5051851862545925, 0.5124998613024105, 0.5199204469598645, 0.5274484767273707, 0.5350855063091983, 0.5428331139349621, 0.5506929006857724, 0.5586664908251056  …  1.7899766970506505, 1.8158941194896647, 1.84218680535361, 1.8688601881549594, 1.8959197800790755, 1.9233711731233287, 1.9512200402527142, 1.9794721365721941, 2.0081333005160213, 2.0372094550542794], iter = 37), [0.001156456913348616 4.7648692493456394e-8 5.456551019037757e-7; 0.00029942249071175333 3.057907732662415e-8 3.49267457951519e-7; … ; 9.407333285062188e-10 0.00016943495830104894 0.0001303664382126568; 9.523762527622712e-10 0.0006498609578498225 0.0005071883069169294], 35)

Just look at the insane speed. Remember we didn’t do any parallelization and stuff. These are just efficient codes.

6.6 Some results from the stationary equilibrium

Share of firms in each state

NX_share, EX_low_share, EX_high_share = sum(pdf_rce, dims = 1)
println("The share of nonexporters is $NX_share, the share of exporters with low technology is $EX_low_share, the share of exporters with high technology is $EX_high_share, in the stationary distribution, ")

The share of nonexporters is 0.4526148278973261, the share of exporters with low technology is 0.3046577047390891, the share of exporters with high technology is 0.24272746736358444, in the stationary distribution, 

plot(ϵ_vec, sol_rce.vNX, xlabel = "productivity", ylabel = "value", label = "Non-exporters")
plot!(ϵ_vec, sol_rce.vEX[:,1], label = "Exporters with low export technology")
plot!(ϵ_vec, sol_rce.vEX[:,2], label = "Exporters with high export technology")

find_cutoff_prod(sol_rce)

3-element Vector{Float64}:
 1.1299695318784664
 0.9239802788519544
 0.7133227740212923

Graphically:

plot(ϵ_vec, sol_rce.polNX, xlabel = "productivity", ylabel = "policy", label = "Non-exporters", yticks = [0,1])
plot!(ϵ_vec, sol_rce.polEX[:,1], label = "Exporters with low export technology")
plot!(ϵ_vec, sol_rce.polEX[:,2], label = "Exporters with high export technology")

plot(ϵ_vec, pdf_rce[:,1], xlabel = "productivity", ylabel = "fraction", label = "Non-exporters", title = "Stationary distribution")
plot!(ϵ_vec, pdf_rce[:,2], label = "Exporters with low export technology")
plot!(ϵ_vec, pdf_rce[:,3], label = "Exporters with high export technology")

Cutoff productivity for firms at different states.

We can collect some moments from the stationary distribution. They will be useful for estimations.

function moments_fixedQ(P, Pstar, sol, pdf, ind_p_mat, ind_pstar_mat; E = E, Qvalue = E.Q, τ = 0.0)
    (;ρϵ, σϵ, grids, ξ) = E

    Π, ϵ_vec = tauchen(grids, ρϵ, σϵ)

    # exporter share
    exporter_share = sum(pdf, dims = 1)[1]

    # starter rate
    starter_rate = sol.polNX⋅pdf[:,1]

    # stopper rate (should be same as starter rate under stationary equilibrium)
    stopper_rate = (1 .- sol.polEX)⋅pdf[:,2:3]
    
    # exporter sales ratio conditional on exporting
    sales_EX_do = zeros(eltype(P), grids, length(ξ))
    sales_EX_fo = similar(sales_EX_do)
    accum = zero(eltype(P))
    for i in 1:grids
        for j in 1:length(ξ)
               sales_EX_do[i, j], sales_EX_fo[i, j]  = sales(1, ϵ_vec[i], j, P, Pstar; E= E, Qvalue = Qvalue, τ = τ)
               accum += sales_EX_fo[i, j]/(sales_EX_do[i, j] + sales_EX_fo[i, j]) * pdf[i, j+1]/(1.0-exporter_share )
        end
    end

    # average export price (in domestic currency) and domestic price ratio conditional on exporting
    fo_price_ratio = sum(Qvalue.*ind_pstar_mat[:,2:3]./ind_p_mat[:,2:3].*pdf_rce[:,2:3])/(1.0-exporter_share)

    return exporter_share, starter_rate, stopper_rate, accum, fo_price_ratio 
end

moments_fixedQ (generic function with 1 method)

moments_fixedQ(P_rce, Pstar_rce, sol_rce, pdf_rce, ind_p_mat_rce, ind_pstar_mat_rce)

(0.19266736559312644, 0.01391198296999171, 0.013911817889600284, 0.46736658942224996, 1.4059666732401315)

Note that by optimality condition of an exporting firm. \(Q p_j^* /p_j = \frac{\xi}{1 - \tau}\), so the last moment, the average export price ratio conditional on exporting can also be computed from the the parameter values of \(\xi\) and \(\tau\) and the conditional fraction of exporters at different states of \(\xi\). This moment can be interpreted as a premium for exporting.

Play with it and try with different parameters. How does the steady state equilibrium change? You might notice something odd when changing the tariff \(\tau\), what drives the result?