Symbolic computation

In [1]:
#Pkg.add("SymPy")
using SymPy, LinearAlgebra
In [2]:
@vars a b c
Out[2]:
(a, b, c)
In [3]:
#we can optionally specify domain with e.g.
x,y,z = symbols("x, y, z", real=true) # also integer, positive, negative...
Out[3]:
(x, y, z)

Algebra

In [4]:
expr = x+2y # note that we do not need to specify 2*y.. this is Julia specific
Out[4]:
\begin{equation*}x + 2 y\end{equation*}
In [5]:
expr2 = x*expr
Out[5]:
\begin{equation*}x \left(x + 2 y\right)\end{equation*}
In [6]:
expanded_expr = SymPy.expand(x*expr)  # Attenction that without the star it gives wrong results!
Out[6]:
\begin{equation*}x^{2} + 2 x y\end{equation*}
In [7]:
factor(expanded_expr)
Out[7]:
\begin{equation*}x \left(x + 2 y\right)\end{equation*}
In [8]:
f = (x+1)^2 # power is ^ and not ** .. again, Julia specific
Out[8]:
\begin{equation*}\left(x + 1\right)^{2}\end{equation*}
In [9]:
g = x^2-2x+1
Out[9]:
\begin{equation*}x^{2} - 2 x + 1\end{equation*}
In [10]:
f-g
Out[10]:
\begin{equation*}- x^{2} + 2 x + \left(x + 1\right)^{2} - 1\end{equation*}
In [11]:
simplify(f-g)
Out[11]:
\begin{equation*}4 x\end{equation*}
In [12]:
expr3 = subs(expr, x=>1, y=>z^3) # substitution can be either numeric or symbolic
Out[12]:
\begin{equation*}2 z^{3} + 1\end{equation*}
In [13]:
solve(expr2,y) # by default solve the equation f(.) = 0
Out[13]:
\[ \left[ \begin{array}{r}- \frac{x}{2}\end{array} \right] \]
In [14]:
solve(expr2,x)
Out[14]:
\[ \left[ \begin{array}{r}0\\- 2 y\end{array} \right] \]
In [15]:
Eq(expr2,2) # create an equation object
Out[15]:
\begin{equation*}x \left(x + 2 y\right) = 2\end{equation*}
In [16]:
solve(Eq(expr2,2),y)
Out[16]:
\[ \left[ \begin{array}{r}- \frac{x}{2} + \frac{1}{x}\end{array} \right] \]
In [17]:
solve([expr,g],[x,y]) # solve systems of equations in multiple variables
Out[17]:
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 [18]:
M = [[1,2,x] [4,5,x]]
Out[18]:
\[\left[ \begin{array}{rr}1&4\\2&5\\x&x\end{array}\right]\]
In [19]:
N = [[x,2] [3,4] [5,6]]
Out[19]:
\[\left[ \begin{array}{rrr}x&3&5\\2&4&6\end{array}\right]\]
In [20]:
M*N
Out[20]:
\[\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 [21]:
det(M*N)
Out[21]:
\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 [22]:
d = [a,b,3]
e = [x,y,2]
Out[22]:
\[ \left[ \begin{array}{r}x\\y\\2\end{array} \right] \]
In [23]:
k = d' * e
Out[23]:
\begin{equation*}x \overline{a} + y \overline{b} + 6\end{equation*}
In [24]:
j = d * e'
Out[24]:
\[\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 [25]:
f = a*x^2-4x+3c
Out[25]:
\begin{equation*}a x^{2} + 3 c - 4 x\end{equation*}
In [26]:
diff(f,x)
Out[26]:
\begin{equation*}2 a x - 4\end{equation*}
In [27]:
solve(diff(f,x),x) # find FOC
Out[27]:
\[ \left[ \begin{array}{r}\frac{2}{a}\end{array} \right] \]
In [28]:
integrate(f,x)
Out[28]:
\begin{equation*}\frac{a x^{3}}{3} + 3 c x - 2 x^{2}\end{equation*}
In [29]:
integrate(f,(x,2,4)) # definite integral between x=2 and x=4
Out[29]:
\begin{equation*}\frac{56 a}{3} + 6 c - 24\end{equation*}

An application: derive Faustmann formula

In [30]:
@vars t k r p NPV       # time, planting costs, interest rate, timber price, Net Present Value (as symbol)
Out[30]:
(t, k, r, p, NPV)
In [31]:
@symfuns S              # generic function for stock
Out[31]:
(S,)

Net Present Value

  • $PV ~=~ \frac{b}{(1+r)^t-1}$
  • $b ~=~ {p S\{t\} - k (1+r)^t}$
In [32]:
NPVf = ( p * S(t) * exp(-r*t)-k)/ (1-exp(-r*t))
Out[32]:
\begin{equation*}\frac{- k + p S{\left (t \right )} e^{- r t}}{1 - e^{- r t}}\end{equation*}
In [33]:
# Derivative of NPV with respect to time
NPVt = diff(NPVf,t)
Out[33]:
\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 [34]:
# Equation of derivative of NPV equal to zero
NPVteq = Eq(NPVt,0)
Out[34]:
\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 [35]:
# 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[35]:
\[ \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 [36]:
# 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[36]:
\begin{equation*}\frac{- k + p S{\left (t \right )} e^{- r t}}{1 - e^{- r t}} = NPV\end{equation*}
In [37]:
# Solve the equation NPV==NPVb with respect to k, which has a unique solution as a function of NPV
k_NPV = solve(NPVeq, k)[1] # k = f(NPV)
Out[37]:
\begin{equation*}\left(- NPV e^{r t} + NPV + p S{\left (t \right )}\right) e^{- r t}\end{equation*}
In [38]:
# Now substitute k
VtNPV = subs(Vt[1], k=>k_NPV) # V' = f(NPV)
Out[38]:
\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 [39]:
# Now simplify
VtNPV2 = simplify(VtNPV) # pS' = f(pS, NPV, r)
Out[39]:
\begin{equation*}r \left(NPV + p S{\left (t \right )}\right)\end{equation*}
In [40]:
# Now expand
SymPy.expand(VtNPV2) # pS' = f(pS, NPV, r), the Faustmann formula
Out[40]:
\begin{equation*}NPV r + p r S{\left (t \right )}\end{equation*}
In [ ]: