\documentclass[../Main.tex]{subfiles}
\graphicspath{{\subfix{Assets/img/}}}

\begin{document}
With the laws of motion introduced in sections \cref{SEC:Laws}, we can now describe
the optimization problem facing each constellation operator.

Each operator recieve utility in each period per 
their per period utility $u^i(\vec s_t,D_t)$, which depends 
on the current sizes of constellations and the level of debris.
In addition, the operator pays for the launch of $x^i_t$ satellites
according to the cost function $F(x)$.
These satellites will become operational in the next period.

Thus the $M$-period (possibly infinite), problem is:
\begin{align}
    \max_{\{\vec x_t\}^M}&~ 
        E\left[ \sum^M_{t=0} \beta^t u^i(\vec s_t, D_t) - F(x^i_t) \right] \\
    &\text{subject to:}\\
    & s^i_{t+1} = (1-l^i(\vec s_t, D_t))s^i_t +x^i_t ~~~  \forall i \\
    & D_{t+1} = (1-\delta)D_t + g(D_t) 
        + \gamma \sum^N_{i=1} l^i(\vec s_t, D_t)
        + \Gamma \sum^N_{i=1} x^i_t
\end{align}
%Assumptions
% - Identical launch costs
% - Identical debris production from destruction.
% 

\subsection{Infinite Period (Bellman) Equation}
The inifinite period version of the problem above can be rewritten in the bellman form as
\begin{align}
    V^i(\vec s_t, \vec x^{\sim i}_t, D_t) = \max_{x^i_t} u^i(\vec s_t, D_t) -F(x) 
    + \beta \left[ V^i(\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1}) \right]
\end{align}
where $x^{\sim i}_t$ represents the launch decisions of all the other constellation
operators.
This implies that the policy function is a best response function, allowing for
a nash equilibrium interpretation of the result.

To solve for the policy function, we have a variety of methods available.
Due to the computational method chosen later, I'm going to examine the conditions
for the existence of an euler equation.


\subsubsection{Euler Equation}
Appendix \cref{APX:Derivations:EulerEquations} contains more details 
on the math involved. 
What follows is just a sketch of the method in matrix notation.

As there is only one choice variable, we get a single optimality condition. 
It can be written in various formats, with the latter matching the appendix the best.
\begin{align}
%    0 =& \parder{}{x^i_t}{} u^i(\vec s_t, D_t) -\parder{}{x^i_t}{}F(x) 
%    + \beta \left[ \parder{}{x^i_t}{} 
%        V^i(\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1}) 
%    \right] \\
    0 =&  -\der{F}{x^i_t}{} 
    + \beta \left[  
        \nabla_{x^i_t} [ \vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1} ]
        \cdot
        \nabla_{\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1}} 
            V^i(\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1}) 
    \right] \label{EQ:OptimalityCondition}\\
    0 =&  -\der{F}{x^i_t}{} + \beta \vec a(\vec s_t,D_t) \cdot \nabla V^i_{t+1} 
    \label{EQ:SimplifiedOptimalityCondition}\\
        =& - f_{x_t} + \beta \vec a_t \cdot \nabla V^i_{t+1}
\end{align}

As there are $N$ constellations we get $N$ satellite stocks, 
$N-1$ decisions $x^{\sim i}$,
and $1$ debris state for a total of $2N$ state 
variables\footnote{recall that $N$ is the number of constellations.}.
Thus there are $2N$ envelope conditions to be found:
\begin{align}
%    \nabla_{\vec s_t, \vec x^{\sim i}_t, D_t} V^i(\vec s_t, \vec x^{\sim i}_t, D_t) 
%    =& \nabla_{\vec s_t, \vec x^{\sim i}_t, D_t} u^i(\vec s_t, D_t) \notag \\
%    &+ \beta \left[ 
%        \nabla_{\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1} }
%            V^i(\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1}) 
%        \cdot
%        \nabla_{\vec s_t, \vec x^{\sim i}_t, D_t}
%            [ \vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1} ]
%    \right] \label{EQ:EnvelopeConditions}
%    \\
    \nabla_{\vec s_t, \vec x^{\sim i}_t, D_t} V^i(\vec s_t, \vec x^{\sim i}_t, D_t) 
    =
    \nabla \vec V^i_t 
    = \vec u^i 
        + \beta B_t \cdot \nabla \vec V^i_{t+1} 
        \label{EQ:SimplifiedEnvelopeConditions}
\end{align}
%When interpreting this, note that 
% $$
% \nabla \vec V^i_{t+1} = \nabla_{[\vec s_{t+1}~ \vec x^{\sim i}_{t+1}~ D_{t+1}] }
%     V^i(\vec s_{t+1}, \vec x^{\sim i}_{t+1}, D_{t+1}) 
% $$
% is a $2N \times 1$ vector of first derivatives but 
% $$
% A = \nabla_{\vec s_t, \vec x^{\sim i}_t, D_t}
%     [ \vec s_{t+1}~ \vec x^{\sim i}_{t+1}~ D_{t+1} ]
% $$
% is a $2N \times 2N$ matrix of first derivatives.

% By solving for $\vec V^i_{t+1}$ as a function of $\vec V^i_{t}$ we get the 
% intertemporal condition:
% \begin{align}
%     \frac{1}{\beta} A^{-1} \left(\nabla \vec V^i_t - \vec u^i_t \right) 
%     = \nabla \vec V^i_{t+1} 
% \end{align}
% Thus one crucial condition for the existence of a solution is that $A^{-1}$ exists for 
% all values the laws of motion and choice functions can take.

% \subsection{Existence}
% I need to do some more diving into conditions for existence.
% Of particular concern is that the way I have specified the debris may lead to
% non-convergence.
% 

To finish constructing the euler equation, we would use the intertemporal 
transition function \cref{EQ:SimplifiedEnvelopeConditions} and iterated
versions of \cref{EQ:OptimalityCondition,EQ:SimplifiedOptimalityCondition}
to construct the $2N+1$ euler equations.\footnote{Double check numbers}
Note that for even a small number of agents -- e.g. 3 -- this iterated substitution
becomes relatively complex, requiring caculating an iterated intertemporal tranisition
function and laws of motion 6 times.
To solve this symbolicly involves inverting a $6 \times 6$ matrix. 
As matrix inversion has approximately an $O(n^3)$ computational complexity,
this becomes unsustainable very quickly.

Section \cref{SEC:Computation} describes how to address this issue to generate
these euler equations using features of modern programming languages and linear algebra
libraries.

\end{document}