Multiple Time Series Analysis

I am reading the book New Introduction to Multiple Time Series Analysis by L├╝tkepohl to refresh my memory of vector autoregressive models, error correction models, and state space models (supposedly all equivalent.) I will build in some computational aspect of these models into this post using python or julia. This post serves as my workbook while reading the book.

To Be Continued…

Make a Dollar, Pave a Road

A popular combinatorics problem is the Making a Dollar problem:

Making a Dollar Given unlimited numbers of 1, 5, 10, 25, and 50-cent denomination coin, how many ways are there to make a dollar.  \blacksquare

The number of ways to mix the coins into a dollar can be counted in two mutually exclusive ways: 1) the number of ways where the biggest denomination is in the coin mix, and 2) the number of ways the biggest denomination is not in the coin mix. If (1) is true, then there is at least one coin of the biggest denomination in the mix. Let us then use one coin of the biggest denomination and count the number of ways to make the remaining amount using using coins of all denominations. If (2) is true, then let us calculate the number of ways to make a dollar using all other denominations. Let f(n,S) be the number of ways to make n-cents using coins with denominations in S = \{d_1, \dots, d_n\}, then

f(n, S) = f(n - d_n, S) + f(n, S \setminus \{d_n\})

This recursive relation needs several boundary conditions: when the amount to make is equal to zero, when the amount is less than zero, and when the set of denominations is empty. The first boundary condition implies that we have just found a way to make the exact change because the remaining amount is zero. The second boundary condition says that we have not been able to found an exact change because the denomination of the coin we just used is too large. The last boundary condition says that we have used up all possible denominations and still cannot make the change. These statements translate respectively to f(0, S) = 1, f(n, S) = 0 for all negative n, and f(n, \varnothing) = 0.

I am picking up a new programming language called Julia so I decided to code this up as an exercise. The procedure described above can be put as follows

function const_comb(sum::Int, parts::Array{Int})
    if sum == 0
        return 1
    end
    if sum < 0 || length(parts) < 1
        return 0
    end
    return const_comb(sum - parts[end], parts) + const_comb(sum, parts[1:end-1])
end

This problem can be generalised in several ways. Suppose it matters what order the coins are returned to make the dollar, for example, we want to treat \{1, 3, 1, 1\} as different from \{1, 1, 1, 3\}, then how many ways can we make a dollar? It may seem necessary to keep track of the sequence to get the changes and then to count each sequence separately. However, this is not too hard to do and we can recycle the algorithm above. Instead of return 1 when we have the the exact change, we return the number of permutations of the coins to get the change, such as the number of permutations of \{1, 1, 1, 3\} to make 6 cents. Such permutation with repeating elements are called a permutation of a multi-set. A multi-set is a collection of objects that need not be distinct, i.e. a multi-set M =\{a_i \cdot d_i, \dots a_k \cdot d_k\} where k, a_i are non-negative integers and d_i are distinct objects, contains a_i of d_i. From combinatorics theory, the number of ways to form a permutation on a multi-set M such that \sum_{i=1}^{k} a_i = a is \frac{a!}{a_1!\dots a_k!}. Sounds like a plan.

Let us first code up the above expression to calculate the number of multi-set permutation. Since it involves a few factorials, special care should be taken to avoid integer overflow. (If you ever wonder why the function lfact–the natural log of the factorial function–needs to exist in Julia’s built-in math library, you will like this post where it shows why a function calculating the hypotenuse of a right triangle need to exist in the cmath library.)

function nPr(n::Int, r::AbstractArray{Int,1})     
    if(any(i->i>n, r))
        return 0
    end
    if(n==0 || all(i->i==0,r))
        return 1
    else
        return exp(lfact(n) - sum(lfact(r)));
    end
end

Now we simply need to keep track of the multi-set of coins of each denomination that make the exact change, and count the number of multi-set permutation. The fact that each object d_i is distinct makes the use of a dictionary in our code a natural choice. The multi-set can be modelled as a dictionary with keys being the distinct object and value being the number of object

function const_comb_rep(sum::Int, parts::AbstractArray{Int, 1}, multiset::Dict)
    if sum == 0
        return nPr(sum(collect(values(multiset))), collect(values(multiset)))
    end
    if sum < 0 || length(parts) < 1
        return 0
    end
    sub_multiset = copy(multiset)
    if haskey(sub_multiset, parts[end])
        sub_multiset[parts[end]] = sub_multiset[parts[end]] + 1
    else
        sub_multiset[parts[end]] = 1
    end
    return const_comb_rep(sum - parts[end], parts, sub_multiset) + const_comb_rep(sum, parts[1:end-1], multiset)
end

running with some test inputs get

const_comb_rep(7, [1, 2, 3], Dict())
44.0
[const_comb_rep(n,1:n, Dict{Int,Int}()) for n = 1:10]
10-element Array{Union{Float64,Int64},1}:
   1.0
   2.0
   4.0
   8.0
  16.0
  32.0
  64.0
 128.0
 256.0
 512.0

As we recurse, the dictionary keeps tract of the number and denomination of the coins used. Finally when we have an exact change, the number of ways to arrange them is computed and returned back.

Note that the number of ways to make Integer n using numbers from 1 \dots n is 2^{n-1}.

Of course, in real life one does not care the order of which changes are returned to make up the right amount. The problem is perhaps more suitable to describe the Pave a Road problem:

Paving a Road Find the number of ways to pave a 1 \times 7 block of road using 1 \times 1, 1 \times 2, and 1 \times 3 blocks, assuming each block is indistinguishable.  \blacksquare

The above algorithm gives 44 ways, and one can actually print out each way if one add a print command inside the if sum == 0 statement.