#Nathan Fox
#Homework 18
#I give permission for this work to be posted online

#Read procedures from class
read(`C19.txt`):

#Help procedure
Help:=proc():
    print(` RandDSP(k, M, U, K) , IsMartingale(P) `):
end:

##PROBLEM 2##

#RandDSP(k, M, U, K): generates a random Discrete-time Stochastic
#process. The inputs is the same as RandMar(k,M,U,K), but the last
#component of the output is not (necessarily) a martingale.
RandDSP:=proc(k, M, U, K) local T, q, NuV, C, G, NonLeaves, i, j, v, Mar, ra:
    T:=RandTree(k, M):
    NuV:=T[1]:
    C:=T[2]:
    G:=T[3]:
    NonLeaves:={seq(op(G[i]), i=1..nops(G)-1)}:
    
    for v in NonLeaves do
        q[v]:=RandLD(nops(C[v+1]), K):
    od:
    
    ra:=rand(1..U):
    for v in G[-1] do
        Mar[v]:=ra():
    od:
    
    for i from k by -1 to 1 do
        for v in G[i] do
            Mar[v]:=ra():
        od:
    od:
    
    return [op(T), [seq(q[i], i=0..nops(NonLeaves)-1)], [seq(Mar[i], i=0..NuV-1)]]:
end:

#IsMartingale(P): inputs a list of length 5 representing a discrete
#stochastic process (of the kind given by the output of procedure
#RandDSP(k, M, U, K)), and outputs true if and only if P is a
#martingale.
IsMartingale:=proc(P) local C, j, L, q, i, a, b:
    C:=P[2]:
    q:=P[4]:
    L:=P[5]:
    
    for i from 0 to nops(C)-1 do
        a:=L[i+1]:
        b:=add(q[i+1][j]*L[C[i+1][j]+1], j=1..nops(C[i+1])):
        if a <> b then
            return false:
        fi:
    od:
    
    return true:
end:

#Test it for random outputs of RandDSP and RandMar (and hopefully get true for the latter).