# Symbolic computation¶

In :
#Pkg.add("SymPy")
using SymPy, LinearAlgebra

In :
@vars a b c

Out:
(a, b, c)
In :
#we can optionally specify domain with e.g.
x,y,z = symbols("x, y, z", real=true) # also integer, positive, negative...

Out:
(x, y, z)

### Algebra¶

In :
expr = x+2y # note that we do not need to specify 2*y.. this is Julia specific

Out:
\begin{equation*}x + 2 y\end{equation*}
In :
expr2 = x*expr

Out:
\begin{equation*}x \left(x + 2 y\right)\end{equation*}
In :
expanded_expr = SymPy.expand(x*expr)  # Attenction that without the star it gives wrong results!

Out:
\begin{equation*}x^{2} + 2 x y\end{equation*}
In :
factor(expanded_expr)

Out:
\begin{equation*}x \left(x + 2 y\right)\end{equation*}
In :
f = (x+1)^2 # power is ^ and not ** .. again, Julia specific

Out:
\begin{equation*}\left(x + 1\right)^{2}\end{equation*}
In :
g = x^2-2x+1

Out:
\begin{equation*}x^{2} - 2 x + 1\end{equation*}
In :
f-g

Out:
\begin{equation*}- x^{2} + 2 x + \left(x + 1\right)^{2} - 1\end{equation*}
In :
simplify(f-g)

Out:
\begin{equation*}4 x\end{equation*}
In :
expr3 = subs(expr, x=>1, y=>z^3) # substitution can be either numeric or symbolic

Out:
\begin{equation*}2 z^{3} + 1\end{equation*}
In :
solve(expr2,y) # by default solve the equation f(.) = 0

Out:
$\left[ \begin{array}{r}- \frac{x}{2}\end{array} \right]$
In :
solve(expr2,x)

Out:
$\left[ \begin{array}{r}0\\- 2 y\end{array} \right]$
In :
Eq(expr2,2) # create an equation object

Out:
\begin{equation*}x \left(x + 2 y\right) = 2\end{equation*}
In :
solve(Eq(expr2,2),y)

Out:
$\left[ \begin{array}{r}- \frac{x}{2} + \frac{1}{x}\end{array} \right]$
In :
solve([expr,g],[x,y]) # solve systems of equations in multiple variables

Out:
1-element Array{Tuple{Sym,Sym},1}:
(1, -1/2)

### Matrix algebra¶

• matrices are just Julian matrices with symbolic entries
• constructing matrices then follows Julia's conventions
• compared with "plain" Julia" here we can do symbolic operations, not just numeric ones
In :
M = [[1,2,x] [4,5,x]]

Out:
$\left[ \begin{array}{rr}1&4\\2&5\\x&x\end{array}\right]$
In :
N = [[x,2] [3,4] [5,6]]

Out:
$\left[ \begin{array}{rrr}x&3&5\\2&4&6\end{array}\right]$
In :
M*N

Out:
$\left[ \begin{array}{rrr}x + 8&19&29\\2 x + 10&26&40\\x^{2} + 2 x&7 x&11 x\end{array}\right]$
In :
det(M*N)

Out:
\begin{equation*}6 x^{2} + 11 x \left(26 x + 208\right) - 11 x \left(38 x + 190\right) - 7 x \left(40 x + 320\right) + 7 x \left(58 x + 290\right) + 12 x\end{equation*}
In :
d = [a,b,3]
e = [x,y,2]

Out:
$\left[ \begin{array}{r}x\\y\\2\end{array} \right]$
In :
k = d' * e

Out:
\begin{equation*}x \overline{a} + y \overline{b} + 6\end{equation*}
In :
j = d * e'

Out:
$\left[ \begin{array}{rrr}a x&a y&2 a\\b x&b y&2 b\\3 x&3 y&6\end{array}\right]$

### Calculus¶

In :
f = a*x^2-4x+3c

Out:
\begin{equation*}a x^{2} + 3 c - 4 x\end{equation*}
In :
diff(f,x)

Out:
\begin{equation*}2 a x - 4\end{equation*}
In :
solve(diff(f,x),x) # find FOC

Out:
$\left[ \begin{array}{r}\frac{2}{a}\end{array} \right]$
In :
integrate(f,x)

Out:
\begin{equation*}\frac{a x^{3}}{3} + 3 c x - 2 x^{2}\end{equation*}
In :
integrate(f,(x,2,4)) # definite integral between x=2 and x=4

Out:
\begin{equation*}\frac{56 a}{3} + 6 c - 24\end{equation*}

## An application: derive Faustmann formula¶

In :
@vars t k r p NPV       # time, planting costs, interest rate, timber price, Net Present Value (as symbol)

Out:
(t, k, r, p, NPV)
In :
@symfuns S              # generic function for stock

Out:
(S,)

### Net Present Value¶

• $PV ~=~ \frac{b}{(1+r)^t-1}$
• $b ~=~ {p S\{t\} - k (1+r)^t}$
In :
NPVf = ( p * S(t) * exp(-r*t)-k)/ (1-exp(-r*t))

Out:
\begin{equation*}\frac{- k + p S{\left (t \right )} e^{- r t}}{1 - e^{- r t}}\end{equation*}
In :
# Derivative of NPV with respect to time
NPVt = diff(NPVf,t)

Out:
\begin{equation*}- \frac{r \left(- k + p S{\left (t \right )} e^{- r t}\right) e^{- r t}}{\left(1 - e^{- r t}\right)^{2}} + \frac{- p r S{\left (t \right )} e^{- r t} + p e^{- r t} \frac{d}{d t} S{\left (t \right )}}{1 - e^{- r t}}\end{equation*}
In :
# Equation of derivative of NPV equal to zero
NPVteq = Eq(NPVt,0)

Out:
\begin{equation*}- \frac{r \left(- k + p S{\left (t \right )} e^{- r t}\right) e^{- r t}}{\left(1 - e^{- r t}\right)^{2}} + \frac{- p r S{\left (t \right )} e^{- r t} + p e^{- r t} \frac{d}{d t} S{\left (t \right )}}{1 - e^{- r t}} = 0\end{equation*}
In :
# Rewriting the equation with the value derivative as dependent variable  (V' = pS')
Vt = solve(NPVteq,p * sympy.Derivative(S(t), t)) # V' = f(pS,r,k,t)

Out:
$\left[ \begin{array}{r}- \frac{r \left(k - p S{\left (t \right )}\right) e^{r t}}{e^{r t} - 1}\end{array} \right]$
In :
# Setting the "symbol" NPV equal to the NPV function
# (because we want to keep in the final output the result in term of NPV)
NPVeq = Eq(NPVf,NPV)

Out:
\begin{equation*}\frac{- k + p S{\left (t \right )} e^{- r t}}{1 - e^{- r t}} = NPV\end{equation*}
In :
# Solve the equation NPV==NPVb with respect to k, which has a unique solution as a function of NPV
k_NPV = solve(NPVeq, k) # k = f(NPV)

Out:
\begin{equation*}\left(- NPV e^{r t} + NPV + p S{\left (t \right )}\right) e^{- r t}\end{equation*}
In :
# Now substitute k
VtNPV = subs(Vt, k=>k_NPV) # V' = f(NPV)

Out:
\begin{equation*}- \frac{r \left(- p S{\left (t \right )} + \left(- NPV e^{r t} + NPV + p S{\left (t \right )}\right) e^{- r t}\right) e^{r t}}{e^{r t} - 1}\end{equation*}
In :
# Now simplify
VtNPV2 = simplify(VtNPV) # pS' = f(pS, NPV, r)

Out:
\begin{equation*}r \left(NPV + p S{\left (t \right )}\right)\end{equation*}
In :
# Now expand
SymPy.expand(VtNPV2) # pS' = f(pS, NPV, r), the Faustmann formula

Out:
\begin{equation*}NPV r + p r S{\left (t \right )}\end{equation*}
In [ ]: