Introduction to Julia

Chase Coleman & Spencer Lyon

3-4-16

Opening Example


In [4]:
# *1                          *2       *3
function bisect(f, a, b, maxit=100, tol::Float64=1e-9)
    fa, fb = f(a), f(b)
    # *4       *5
    for it in 1:maxit 
        mid = (a + b)/2
        fmid = f(mid)
        
        # *6
        if abs(fmid) < tol            
            # *7
            return mid
        end

        if fa*fmid > 0
            fa, a = fmid, mid  
        else
            fb, b = fmid, mid
        end
    end
    
    # *8
    error("maximum iterations exceeded")
end


Out[4]:
bisect (generic function with 3 methods)
  1. Define new functions with function ... end
  2. Default arugments (..., arg=default_value)
  3. Typed arguments (..., arg::Type)
  4. For loop for X (in|=) SOMETHING ... end
  5. Create ranges A:B (not dense like Matlab)
  6. If statement if CONDITION BLOCK end
  7. Return statement return STUFF (optional, see next example)
  8. Throwing error error(MESSAGE)

In [5]:
# shorthand function syntax
f(x) = x^2-2

# longer syntax --  equivalent to above
function f2(x)
    x^2-2
end

# much longer syntax -- still equivalent to above
function f3(x)
    return x^2-2
end


Out[5]:
f3 (generic function with 1 method)

In [6]:
println(bisect(f, 0, 4))
println(bisect(f2, 0, 4))
println(bisect(f3, 0, 4))


1.4142135623842478
1.4142135623842478
1.4142135623842478

Types

Everything in Julia has a type.

You can inspect the type of an object using the typeof function:


In [18]:
for obj in [1.0, 1, "foo", "φ", bisect, Int8(4), Float64, true]
    # Notice string interpolation syntax `$`
    println("$obj is a $(typeof(obj))")
end


1.0 is a Float64
1 is a Int64
foo is a ASCIIString
φ is a UTF8String
bisect is a Function
4 is a Int8
Float64 is a DataType
true is a Bool

Type Parameters

Types can have parameters.

This concept is best understood by example


In [20]:
x = [1, 2, 3]
typeof(x)  # whats the `{` and `}` stuff all about?


Out[20]:
Array{Int64,1}

Int and 1 are type parameters. In this case it tells us that the array is filled with Ints and has 1 dimension (it is a vector, yes Julia has vectors).

Allow you to do all sorts of magic that will become clear later. For now just recognize the { } syntax.

User defined types

You can define your own types (you should do this a lot).

Two forms of types:

  • abstract: you can't create these, but they help you group related types together
  • composite: you do create these, they are the actual data of your program

In [21]:
abstract Exog

type AR1 <: Exog
    rho::Float64  # you should put types on the fields of your types
    sigma::Float64
end

"""
I'm docstring. I describe the `MarkovChain` type.

`x0` is the initial distribution
"""
type MarkovChain{T} <: Exog
    Π::Matrix{Float64}
    vals::AbstractVector{T}
    x0::Vector{Float64}
end

# functions to give `vals` and `x0` default arguments
MarkovChain(Π, v) = MarkovChain(Π, v, fill(1/length(v), length(v)))
MarkovChain(Π) = MarkovChain(Π, 1:size(Π, 1))

In [24]:
AR1(0.9, 0.1)


Out[24]:
AR1(0.9,0.1)

In [23]:
?MarkovChain


search: 
Out[23]:

I'm docstring. I describe the MarkovChain type.

x0 is the initial distribution

MarkovChain

Multiple Dispatch

Function behavior can be specialized based on the type (and number) of all function arguments

Let's see some examples


In [41]:
g(x) = "I have something"
g(x::Int) = "I have an integer"
g(x::Float64) = "I have a float"
g(x::Number) = "I have some kind of number"
g(x::Array) = "I have an array"

for x in ("hello", 1, 1.0, 1//2, [1, 2, 3])
    @printf "%-12s%s\n" "g($x)" g(x)
end


g(hello)    I have something
g(1)        I have an integer
g(1.0)      I have a float
g(1//2)     I have some kind of number
g([1,2,3])  I have an array

We can add methods to the g function that take multiple arguments

Notice how the return value depends on the types of both arguments


In [51]:
g(x, y) = "I have two things"
g(x::Int, y) = "I have an integer and something else"
g(x::Int, y::Number) = "I have an integer and a number"
g(x::Int, y::Int) = "I have two integers"
g(x::Array, y::Array) = "I have two arrays"
g(x::Array{Float64}, y::Array{Float64}) = "I have two arrays that have floats"


Out[51]:
g (generic function with 13 methods)

In [52]:
stuff = (("x", "y"), (1, "x"),  (1, 1//2), 
         (1, 2),  (1, 2.0), ([1], [2]), 
         ([1.0], [2.0]))

for (x1, x2) in stuff
    @printf "%-18s%s\n" "g($x1, $x2)" g(x1, x2)
end


g(x, y)           I have two things
g(1, x)           I have an integer and something else
g(1, 1//2)        I have an integer and a number
g(1, 2)           I have two integers
g(1, 2.0)         I have an integer and a number
g([1], [2])       I have two arrays
g([1.0], [2.0])   I have two arrays that have floats

Example

Let's construct routines that will allow us to simulate any subtype of Exog

To do this we will need each subtype of Exog to implement a iter method

This method should take two arguments:

  • The Exog subtype
  • The current state

It should return the state on the next iteration


In [58]:
iter(ar1::AR1, x) = ar1.rho*x + ar1.sigma*randn()

function iter(mc::MarkovChain, s::Int, v)
    ind = searchsortedfirst(cumsum(vec(mc.Π[s, :])), rand())
    return mc.vals[ind]
end

iter{T}(mc::MarkovChain{T}, v::T) = iter(mc, findfirst(mc.vals, v), v)
iter(::Exog, x) = error("iter should be implemented by each Exog subtype")


Out[58]:
iter (generic function with 4 methods)

Now we will define single simulate function for all Exog subtypes


In [64]:
# NOTE `;` for keyword argument
function simulate(ex::Exog, x0; capT::Int=10)
    out = Array(typeof(x0), capT)
    out[1] = x0

    for t = 2:capT
        out[t] = iter(ex, out[t-1])
    end
    out
end

# for MarkovChain we have more info, so we don't need to give it an x0
# define another method that hands off to the method above
function simulate(mc::MarkovChain; capT::Int=100)
    v = mc.vals[searchsortedfirst(cumsum(mc.x0), rand())]
    simulate(mc, v; capT=capT)
end


Out[64]:
simulate (generic function with 2 methods)

In [70]:
ar1 = AR1(0.9, 0.1)
simulate(ar1, 0.5)


Out[70]:
10-element Array{Float64,1}:
 0.5     
 0.415777
 0.364054
 0.320674
 0.300484
 0.456054
 0.316475
 0.303914
 0.35384 
 0.282838

In [69]:
mc = MarkovChain([0.7 0.3; 0.4 0.6], [0.5, 0.9], [0.2, 0.8])
simulate(mc; capT=5)


Out[69]:
5-element Array{Float64,1}:
 0.5
 0.5
 0.5
 0.9
 0.5