Intro

Trying stuff



In [2]:

    
using Plots
ts_length = 100
epsilon_values = rand(ts_length)
plot(epsilon_values, color="blue")









    Out[2]:

The effect of the statement using Plots is to make all the names exported by the Plots module available in the global scope.

If you prefer to be more selective you can replace using Plosts with: import Plots: plot

Now only the plot function is accessible.

Since our program uses only the plot function from this module, either would have worked in the previous example.

Arrays

The function call epsilon_values = randn(ts_length) creates one of the most fundamental Julia data types: an array.



In [3]:

    
typeof(epsilon_values)









    Out[3]:





Array{Float64,1}



In [4]:

    
epsilon_values









    Out[4]:





100-element Array{Float64,1}:
 0.134494 
 0.792362 
 0.817166 
 0.570655 
 0.849157 
 0.254228 
 0.296585 
 0.867184 
 0.935106 
 0.315811 
 0.0853013
 0.511376 
 0.591898 
 ⋮        
 0.0912906
 0.402971 
 0.176126 
 0.616558 
 0.420713 
 0.586656 
 0.905527 
 0.229745 
 0.296782 
 0.354149 
 0.411254 
 0.620829

The information from typeof() tells us that epsiolon_values is an array of 64 bit floating point values, of dimension 1.

Julia arrays are quite flexible — they can store heterogeneous data for example:



In [5]:

    
x = [10, "foo", false]









    Out[5]:





3-element Array{Any,1}:
    10     
      "foo"
 false

Notice now that the data type is recorded as Any, since the array contains mixed data. The first element of x is an integer.



In [6]:

    
typeof(x[1])









    Out[6]:





Int64

The second is a string.



In [7]:

    
typeof(x[2])









    Out[7]:





String

The third is the boolean value false.



In [8]:

    
typeof(x[3])









    Out[8]:





Bool

Notice from the above that

array indices start at 1 (unlike Python, where arrays are zero-based)
array elements are referenced using square brackets (unlike MATLAB and Fortran)

Julia contains many functions for acting on arrays — we'll review them later.

For now here's several examples, applied to the same list x = [10, "foo", false]



In [9]:

    
length(x)









    Out[9]:





3



In [10]:

    
pop!(x)









    Out[10]:





false



In [11]:

    
x









    Out[11]:





2-element Array{Any,1}:
 10     
   "foo"



In [12]:

    
push!(x, "bar")









    Out[12]:





3-element Array{Any,1}:
 10     
   "foo"
   "bar"



In [13]:

    
x









    Out[13]:





3-element Array{Any,1}:
 10     
   "foo"
   "bar"

The first example just returns the length of the list.

The second pop!(), pops the last element off the list and returns it.

In doing so it changes the list (by dropping the last element).

Because of this we call pop! a mutating method.

It's conventional in Julia that mutating methods end in ! to remind the user that the function has other effects beyond just returning a value.

The function push!() is similar, except that it appends its second argument to the array.

For Loops

Although there's no need in terms of what we wanted to achieve with our program, for the sake of learning syntax let's rewrite our program to use a for loop.



In [14]:

    
ts_length = 100
epsilon_values = Array{Float64}(ts_length)
for i in 1:ts_length
    epsilon_values[i] = randn()
end
plot(epsilon_values, color="red")









    Out[14]:

Here we first declared epsilon_values to be an empty array for storing 64 bit floating point numbers.

The for loop then populates this array by successive calls to randn().

Called without an argument, randn() returns a single float

Like all code blocks in Julia, the end of the for loop code block (which is just one line here) is indicated by the keyword end.

The word in from the for loop can be replaced by symbol =.

The expression 1:ts_length creates an iterator that is looped over — in this case the integers from 1 to ts_length.

Iterators are memory efficinet because the elements are generated on the fly rather than stored in memory.

In Julia you can also loop directly over arrays themselves, like so



In [15]:

    
words = ["foo", "bar"]
for word in words
    println("Hello $word")
end









    



Hello foo
Hello bar

While Loops

The syntax for the while loop contains no surprises



In [16]:

    
ts_length = 100
epsilon_values = Array{Float64}(ts_length)
i = 1
while i <= ts_length
    epsilon_values[i] = randn()
    i = i + 1
end
plot(epsilon_values, color="green")









    Out[16]:

The next example does the same thing with a condition and the break statement.



In [17]:

    
ts_length = 100
epsilon_values = Array{Float64}(ts_length)
i = 1
while true
    epsilon_values[i] = randn()
    i = i + 1
    if i > ts_length
        break
    end
end
plot(epsilon_values, color="black")









    Out[17]:

User-Defined Functions

For the sake of the exercise, let's now go back to the for loop but restructure our program so that generation of random variables takes place within a user-defined function.



In [18]:

    
function generate_data(n)
    epsilon_values = Array{Float64}(n)
    for i in 1:n
        epsilon_values[i] = randn()
    end
    return epsilon_values
end

ts_length = 100
data = generate_data(ts_length)
plot(data, color="brown")









    Out[18]:

Here

function is a Julia keyword that indicates the start of a function definition
generate_data is an arbitrary name for the function
return is a keyword indicating the return value

A Slightly More Useful Function

Of course the function generate_data is completely contrived.

We could just write the following and be done:



In [19]:

    
ts_length = 100
data = randn(ts_length)
plot(data, color="grey")









    Out[19]:

Let's make slightly more useful function.

This function will be passed a choice of probability distribution and respond by plotting a histogram of observations.

In doing so we'll make use of the Distributions package



In [20]:

    
Pkg.add("Distributions")









    



INFO: Package Distributions is already installed
INFO: METADATA is out-of-date — you may not have the latest version of Distributions
INFO: Use `Pkg.update()` to get the latest versions of your packages

Here's the code



In [21]:

    
using Distributions

function plot_histogram(distribution, n)
    epsilon_values = rand(distribution, n)  # n draws from distribution
    histogram(epsilon_values)
end

lp = Laplace()
plot_histogram(lp, 500)









    Out[21]:

Let's have a casual discussion of how all this works while leaving technical details for later in the lectures.

First, lp = Laplace() creates an instance of a data type defined in the Distributions module that represents the Laplace distribution.

The name lp is bound to this object.

When we make the function call plot_histogram(lp, 500) the code in the body of the function plot_histogram is run with

the name distribution bound to the name object as lp
the name n bound to the integer 500

A Mystery

Now consider the function call rand(distribution, n)

This looks like something of a mystery.

The function rand() is defined in the base library such that rand(n) returns n uniform random variables on [0,1)



In [22]:

    
rand(3)









    Out[22]:





3-element Array{Float64,1}:
 0.905887 
 0.101278 
 0.0489148

On the other hand, distribution points to a data type representing the Laplace distribution that has been defined in a third party package.

So how can it be that rand() is able to take this kind of object as an argument and return the output that we want?

The answer in a nutshell is multiple dispatch

This refers to the idea that functions in Julia can have different behaviour depending on the particular arguments that they're passed.

Hence in Julia we can take an existing function and give it a new behaviour by defining how it acts on a new type of object.

The interpreter knows which function definition to apply in a given setting by looking at the types of the objects the function is called on.

In Julia these alternative versions of a function are called methods

Exercises

Exercise 1

Recall that $n!$ is read as "n factorial" and defined as $n! = n×(n-1)×…×2×1$

In Julia you can compute this value with factorial(n).

Write your own version of this function called factorial2, using for loop.



In [33]:

    
function factorial2(n)
    k = 1
    for i in 1:n
        k = k * i
    end
    return k
end

print(factorial(4),'\n',factorial2(4))

Exercise 2

The binomial random variable $Y \sim Bin(n,p)$ represents:

number of successes in $n$ binary trials
each trial succeeds with probability $p$

Using only rand() from the set of Julia's build-in random number generators (not the Distributions package), write a function binomial_rv such that binomial_rv(n, p) generates one draw of $Y$

Hint: If $U$ is uniform on (0,1) and $P \in (0,1)$, then the expression U < p evaluates to true with probability p



In [35]:

    
function binomial_rv(n, p)
    count = 0
    U = rand(n)
    for i in 1:n
        if U[i] < p
            count += 1
        end
    end
    return count
end

for j in 1:25
    b = binomial_rv(10, 0.5)
    print("$b, ")
end









    



5, 3, 6, 5, 4, 3, 9, 3, 4, 6, 6, 5, 5, 4, 2, 6, 4, 4, 6, 5, 3, 6, 5, 4, 6,

Exercise 3

Compute an approximation to $\pi$ using Monte Carlo.

For random number generations use only rand()

Your hints are as follows:

If $U$ is a bivariate uniform random variable on the unit square $(0,1)^2$, then the probability that $U$ lies in a subset $B$ of $(0,1)^2$ is equal to the area of $B$
If $U_1,...,U_n$ are iid copies of $U$, then, as $n$ gets large, the fraction that falls in $B$ converges to the probability of landing in $B$
For a circle, $A = \pi r^2$



In [38]:

    
n = 10000000

count = 0
for i in 1:n
    u, v = rand(2)
    d = sqrt((u - 0.5)^2 + (v - 0.5)^2)  # Distance from middle of square
    if d < 0.5
        count += 1
    end
end

area_estimate = count / n

print(area_estimate * 4)  # Dividing by radius**2

Exercise 4

Write a program that prints one realization of the following random device:

Flip an unbiased coin 10 times
If 3 consecutive heads occur one or more times within this sequence, pay one dollar
If not, pay nothing

Once again use only rand() as your random number generator



In [57]:

    
payoff = 0
count = 0

print("Count = ")

for i in 1:10
    x = rand()
    if x > 0.5
        count += 1
    else
        count = 0
    end
    print(count)
    
    if count >= 3
        payoff += 1
    end
end


print("\nPayoff = $payoff")









    



Count = 0012300120
Payoff = 1



In [58]:

    
payoff = 0
count = 0

print("Count = ")

for i in 1:10
    x = rand()
    count = x < 0.5 ? count + 1 : 0
    print(count)
    
    if count >= 3
        payoff += 1
    end
end

print("\nPayoff = $payoff")









    



Count = 1234012300
Payoff = 3

Exercise 5

Simulate and plot the correlated time series.

$$x_{t+1} = \alpha x_t + \varepsilon_{t+1} \quad \textrm{where} \quad x_0 = 0 \quad \textrm{and} \quad t = 0,...,T$$

The sequence of shocks $\{\varepsilon_t\}$ is assumed to be iid and standard normal

Set $T = 200$ and $\alpha = 0.9$



In [63]:

    
# Set parameters
alpha = 0.9
T = 200

# Make an array of 200 zeroes
x = zeros(T + 1)

# Construct `for` loop
for t in 1:T
    x[t+1] = alpha * x[t] + randn()
end

# Plot graph
plot(x)









    Out[63]:

Exercise 6

Plot three simulated time series, one for each of the cases $\alpha = 0$, $\alpha = 0.8$ and $\alpha = 0.98$

In particular, you should produce (modulo randomness) a figure that looks as follows:

(The figure illustrates how time series with the same one-step-ahead conditional volatilities, as these three processes havem can have very different unconditional volatilities)



In [67]:

    
alphas = [0.0, 0.8, 0.98]
T = 200

series = []
labels = []

for alpha in alphas
    x = zeros(T + 1)
    x[1] = 0
    for t in 1:T
        x[t+1] = alpha * x[t] + randn()
    end
    push!(series, x)
    push!(labels, "alpha = $alpha")
end

plot(series, label=reshape(labels, 1, length(labels)))









    Out[67]:



In [68]:

    
?reshape









    



search: reshape promote_shape InverseWishart







    Out[68]:





reshape(A, dims...) -> R
reshape(A, dims) -> R
Return an array R with the same data as A, but with different dimension sizes or number of dimensions. The two arrays share the same underlying data, so that setting elements of R alters the values of A and vice versa.
The new dimensions may be specified either as a list of arguments or as a shape tuple. At most one dimension may be specified with a :, in which case its length is computed such that its product with all the specified dimensions is equal to the length of the original array A. The total number of elements must not change.

jldoctest
julia> A = collect(1:16)
16-element Array{Int64,1}:
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16

julia> reshape(A, (4, 4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> reshape(A, 2, :)
2×8 Array{Int64,2}:
 1  3  5  7   9  11  13  15
 2  4  6  8  10  12  14  16



In [ ]: