Last Update: May 18, 2015.
We will first learn Maple, and how to program with it. This semester we will learn how to approach Probability and Statistics the way it SHOULD have been developed, had the computer been around at the time of Pascal and Fermat, using what I call "symbol-crunching". After mastering the foundation of probability, we will learn Financial Calculus, that happened to be a very beautiful part of mathematics (even if you are not fond of money), using our symbolic-computational approach.
In addition to the actual, very important content, students will master the methodology of computer-generated and computer-assisted research that is so crucial for their future.
There are no prerequisites (in particular, no knowledge of probability, Maple, or programming experience is assumed). Also, no overlap with previous years.
Make sure that you are in your home directory, and type the following
Subject: ExpMathRocks; My homework directory
with the url of your homework directory. For example,
http://sites.math.rutgers.edu/~nhf12/em15
[For experts only] Generalize DoN(N) and write a program
DoNg(N,p)
where the probability of winning a dollar is p.
You should assume that p is a rational number.
#Name
#Date, Assignment hw2.txt
#OK to Post
OR:
#Please Do Not Post
Hint: to generate a loaded coin with probability p=a/b, where and b are integers,
use
evalb(rand(1..b)()<=a);
Hint: instead of the set of birthdays, consider the list of birthdays, and then sort them [using the Maple command sort] and check if no consecutive entries are ≤ r apart]
Check your formula (approximately, of course), by running the simulation program BDg(k,N,r) many times for various k,N,r. In particular for k=19, N=365, and r=1.
IsDerangement(pi)
that inputs a permutation pi, and outputs true if the permutation is a derangement (i.e. has no fixed points) and false otherwise. For example IsDerangement([2,3,4,1]) should be true but IsDerangement([2,4,3,1]) should be false.
contains Procedures
MS(N,k,M): running S(N,k) M times and recording the smallest estimate and the largest
Contains procedures
Pr(B|A)= Pr(A|B)Pr(B)/Pr(A)
#Name
#Date, Assignment hw3.txt
#OK to Post
OR:
#Please Do Not Post
MS(N,k,M,eps)
That inputs a positive integer N, a positive integer k (smaller than N, in practice much smaller), a large integer M, and a small number eps, and outputs the number of times the estimated probability was not in the interval [1/2-eps,1/2+eps]
What did you get for
PostV(L,M)
that checks that L and M are good inputs, and returns FAIL if it is not, and prints out an explanation why it is not. Remember that L is a list, M is a list of lists each of its elements has the same length, L can only have probabilities, they should add up to 1, etc. etc.
contains Procedures:
#Name
#Date, Assignment hw4.txt
#OK to Post
OR:
#Please Do Not Post
ExtractData(A,i,j,a)
that inputs a list of lists of statistical data, like in ClassStat() of C4.txt and outputs the sublist consisting only of those items for which the i-th component equals a, and lists only the j-th component. For example
ExtractData([[John,M,25,180,Bl],[Jane,F,30,170,Br],[Jim,M,34,190,Bl],[Nancy,F,50,165,Bl]],2,4,M);
should output
[180,190]
Use this procedure, as well as procedures AveS(L),and sd(P,L) (that we did in class) to find
combinat[permute](n)
(for the desired n), for n from 1 to 7. Can you conjecture an explicit expression for these from the seven data points?
[Hint: first write a procedure, fp(pi), that inputs a permutation pi and outputs its number of fixed points
(i.e. the number of i's such that pi[i]=i),for example fp([3,2,1,4])=2.
HowManyGames(A,m,n)
that inputs a set of positive integers, A, and two positive integers, m and n, and outputs the number of possible scoring histories that could have lead to the score m:n in a game where the `atomic' scores are taken from the set A. For example,
HowManyGames({1},4,3)
tells you the number of ways a Soccer game whose final score was 4:3 could have evolved.
[Hint: the number of possible game histories with score m:n is the coefficient of x^m*y^n in
the bi-variate Maclaurin expansion of
1/(1-add(x^a, a in A)-add(y^a, a in A))
]
Consider the enumerating sequence for games that ended in a tie
seq(HowManyGames(A,i,i),i=0..10)
For which games is the sequence already in the OEIS. You are welcome to submit those that are not
contains Procedures (in addition to those of C4.txt):
#Name
#Date, Assignment hw5.txt
#OK to Post
OR:
#Please Do Not Post
Hint: the probability generating function is (p*x+(1-p)/x)^n . Another hint, since square-roots are a pain, it is more convenient to consider, for the odd standardized moments, their square, take the limit of the square, and then take the square-root (the sqrt of 0 is 0).
SymDie(f,n,k)
that inputs a specific positive integer, f, a symbol (or number) n, and a positive integer k larger than 2, and outputs the expectation, variance, and third through k-th moment of the random variable "total number of dots" if you roll a FAIR f-faced die n times. Using the output for various f,
CP(P1,L1,P2,L2)
that inputs
1. a finite probability distribution P1, of length n1, say, and a r.v. L1, of length n1,
2. a finite probability distribution P1, of length n2, say, and a r.v. L2, of length n2,
and outputs the pair [P,L] where P is a list of size n1*n2 where the probability distribution defined on pairs [i,j] (1 ≤ i ≤ n1, (1 ≤ j ≤ n2), (order this set of pairs in lexicographic order), and L is the r.v. defined by
L[i,j]:=L1[i]+L2[j]
For example
CP([1/2,1/3],[5,8],[1/10,9/10],[-2,3]) ;
should output (with the ordering of pairs [[1,1],[1,2],[2,1],[2,2]] [previous version had an error, spotted by Melissa Romanus]
[1/20,9/20,1/30,9/30],[3,8,6,11]]
Experiment with various random choices of inputs, and verify, experimentally that
PGF(op(CP(P1,L1,P2,L2)),x)=PGF(P1,L1,x)*PGF(P2,L2,x)
a1,a2,a3,a4,a5,a6 (all POSITIVE inetegers, repeats are allowed)
and
b1,b2,b3,b4,b5,b6 (all POSITIVE inetegers, repeats are allowed)
such that if you throw them, you would have the same probability of getting the sum to be k as if you threw the usual 2 dice, for every k from 2 to 12. In other words such that
PGF(op(CP([(1/6)$6],[a1,a2,a3,a4,a5,a6],[(1/6)$6],[b1,b2,b3,b4,b5,b6])),x)= PGF(op(CP([(1/6)$6],[1,2,3,4,5,6],[(1/6)$6],[1,2,3,4,5,6])),x)
contains Procedures
#Name
#Date, Assignment hw6.txt
#OK to Post
OR:
#Please Do Not Post
[Corrected Feb. 12, 2015, thanks to Aniket Shah]
The classical significance tests use the Central Limit Theorem that uses the normal distribution,
exp(-x^2/2)/sqrt(2*Pi) , -infinity < x < infinity
and the error function
Φ(x)=int(exp(-t^2/2)/sqrt(2*Pi), t=-infinity..x)
as follows. Given the random binomial distribution, B(n,p) with expectation np and variance np(1-p), form the standardized r.v.,
Y:=(X-np)/sqrt(n*p*(1-p))
Pr(a ≤ Y ≤ b) is approximated, thanks to CLT, by
Φ(b)-Φ(a)=int(exp(-t^2/2)/sqrt(2*Pi), t=a..b)
Write a procedure
AcceptNHnormal(p,n,k,eps)
that does the same thing as AcceptNH(p,n,k,eps), but accepts the Null Hypothesis if
Φ(-abs(Y)) ≥ eps/2
Experiment with several n,k, and p, and find out whether you get the same output.
Explain why you get the same thing (or almost).
By applying ID(P^n,x,k), with k=6, and f a GENERAL (SYMBOLIC!) pgf
P=p[1]*x^(L[1])+p[2]*x^(L[2])+p[3]*x^(L[3])+p[4]*x^(L[4]) ,
and using the Maple command
lim( ...., n=infinity)
Prove through the 6-th moment, the ASYMPTOTIC NORMALITY, of the random variable "total sum of dots" obtained by rolling n times, a FOUR-faced (i.e. tetrahedral), not necessarily fair, die.
Pr(P,x,n,a,b)
that inputs an integer n, and integers a and b, computes the probability that if the die whose probability generating function (for one roll) is the polynomial P(x) in x, is rolled n times, the total sum of the dots is between a and b
Hint: Expand P(x)^n and collect (using the add command) all the coefficients of x between x^a and x^b
Pr(P,x,n, E+c1*sqrt(V),E+c2*sqrt(V)) and Φ(c2)-Φ(c1)
[Using the E and V that you found above]
for
P=(x+2*x^2+x^3)/4 c1=-2 c2=2 and for n=200
contains Procedures
(x^2+y^2+2*a*y)^2=2*a^2*(x^2+y^2)
( a is a parameter, experiment with a=1, a=2, etc. and see what looks best)
WITHOUT axes, (axes=none), also use the textblock command to insert, inside this curve, the phrase
YourName LOVES EXPERIMENTAL MATHEMATICS
where YourName is replaced by your name. If possible, have colors.
By using help, (or Garvan's book), find out how to save your plot as .jpg or .gif file.
#Name
#Date, Assignment hw7.txt
#OK to Post
OR:
#Please Do Not Post
BPg(p,n,k,c0,a,b)
where c0 is the prior probability that the creator of the coin actually made the probability of a Head p, and if not, with prob. 1-c0, he chooses, uniformly at random, a probability from the interval [a,b]. For example BPg(p,n,k,1/2,0,1) should give the same output as BP(p,n,k) in file C7.txt .
Experiment with various reasonable inputs, for example
BPg(1/2,1000,550,2/3,.4,.6) , BPg(1/3,3000,1100,1/4,.3,.7), etc.
Write a procedure
OneRun(p,n)
that simulates tossing, n times, a loaded coin whose probability of Heads in p (assume that p is a rational number a/b, so that you can use the rand(1..b) command to simulate a coin toss with Prob. a/b), and returns the number of times it got Heads.
ENHnS(p,n,eps,K): that runs the above-mentioned procedure OneRun(p,n) K times (K at least 100), looks at its output, and, using ENHn(p,n,eps), decides whether to accept the Null Hypothesis, and record the ratio of those where you accepted it (correctly, of course).
Try out ENHnS(1/3,300,1/20,1000)
several times, and see whether you get numbers close to .95
SimulatePoisson(n,a,K)
where n is a positive integer at least 10, and K a positive integer at least 100, and a is a small integer, runs the procedure
OneRun(a/n,n)
K times
and outputs the list L=[L0,L1,L2,L3], where L0 is the ratio (out of the K trials) that yielded 0, L1 is the ratio (ouf ot the trials that yielded 1), etc. Does this agree with the theoretical values?
Hint (for novices): define a table T[i], for i=0,1,2,3, initialize them all to 0, and then whenever you get an output from OneRun(a/n,n), call it k, and update T[k] to be T[k]+1. When you are done with the K runs, summarize the output, in floating point as [T[0]/K,T[1]/K,T[2]/K,T[3]/K]]
contains Procedures
#Name
#Date, Assignment hw8.txt
#OK to Post
OR:
#Please Do Not Post
f(N,a):=the probability of exiting a winner (with N dollars) if you start out with a dollars, and at every round you win a dollar with probability 1/2 and you lose a dollar with probability 1/2.
Having conjectured that expression, prove it rigorously [Hint: the system of linear equations in procedure PrE uniquely determines the probabilities, so once you have a conjectured expression, all you have to do is plug-it in and check the two boundary values.]
Modify procedure PrE(N,p), and write a procedure
ExpectedDurationE(N,p)
that inputs N (the exit capital if you are a winner), and p (the probability of winning a dollar at each round), and outputs the list of length N-1 whose a-th entry is the expected duration of the game (exiting either a winner or loser).
[Hint: Set up a system of linear equations as before, but remember that if a=0 or a=N then the game is over so the boundary conditions are g(0)=0 g(N)=0. If you currently have a dollars, then with probability p you have a+1 next round, and probability 1-p, you have a-1 dollars next round, but don't forget to account for the current round, your remaining number of days, today, left to you until you die,equals one more than the number of days left to you, tomorrow. ]
[Added Feb. 17, 2015, 9:15am: I thank John Kim for spotting a typo in the previous version]
ExpectedDurationE(N,1/2)
guess an explicit formula for the expected duration of a game with maximal capital N, starting with a dollars, in the fair case. Let's call it g(a,N). Having guessed it, prove it rigorously!
GRsim(N,a,1/2,1000)[2]
for various N and a.
PrE(N,9/19)[N/2]
for N=50,80,100
[Added Feb. 17, 2015, 9:15am: I thank John Kim and Melissa Romanus for realizing that the previous values of N=1000, N=2000, were too big for our modest computers]
GRsim(N,a,p,K):
and call it
GRsimWL(N,a,p,K):
that returns a triple of numbers. The first number is as before, the ratio of wins to K, the second number is the average duration of the games where he exited a winner, and the third number is the average duration of the games where he exited a loser.
contains Procedures (in addition to the ones from C8.txt)
#Name
#Date, Assignment hw9.txt
#OK to Post
OR:
#Please Do Not Post
[Hint, it is more convenient to first write a more general procedure that inputs an arbitrary list L of positive integers and outputs the set of lists of positive integers that are pointwise smaller-or-equal than L. You may use a recursive algorithm, that chops the last entry of L, finds the corresponding smaller set of shorter lists, and then extends the lists]
RankStrategies(p,N,i)
That inputs a probability p, a positive integer N (denoting, as usual the exit capital), and an integer i between 1 and N-1, denoting the starting capital, and outputs the ranking of the strategies (arranged as a list of lists) from best to worst, if you play the gambling game with probability p, and maximal capital N.
For N=8, and p=9/19, do you get the same output (ranking) for all initial capitals i from 1 to 7? How about N=8 and p=10/19?
RankStrategies6(p,i)
That inputs a probability p, and an integer i between 1 and 5, denoting the initial capital when you entered the casino, and outputs the ranking of the strategies (arranged as a list of 1*2*3*2*1=12 lists of length 5) from best to worst if you play the gambling game with probability p, and maximal capital 5.
Do you get the same output (ranking) for all initial capitals i from 1 to 5, where p=9/19? Where p=10/19?
ExpectedDurationE(N,p)
of Homework assignment 8 and write a procedure
ExpectedDurationEG(N,p,L)
that handles, a game with a general strategy L.
Contains Procedures:
#Name
#Date, Assignment hw10.txt
#OK to Post
OR:
#Please Do Not Post
PGL(p,N,x)
that inputs a numeric or symbolic probability p, a numeric positive integer N, and a variable x, and outputs the list of length N-1 whose n-th entry is the probability generating function (that must be a rational function of x (why?)) for the random variable: `duration of game' started with n dollars, ending with 0 or N, and with probability of a win in each round p.
By using simplify and/or normal, make sure that PGL(p,N,x) agrees with the output of PGFL(p,n,N,x) for N=10 and N=16.
ID(f,x,k)
from C6.txt , and write a procedure
IDlimit(f,x,k)
IDlimit(PGFL(1/2,n,2*n,x),x,4);
[Warning: be patient, it may take a few minutes], and rederive RIGOROUSLY the fact that
E[Xn]=n2
More impressively, derive (or rather let your computer derive) an explicit (polynomial) expression in n, for the Variance, and convince yourself that there is no concentration about the mean. Even more impressively, find an explicit expression in n, for the standardized third and fourth moments, and find their limits as n goes to infinity. Is Xn asympotically normal, as n goes to infinity?
Contains (only one) procedure:
#Name
#Date, Assignment hw11.txt
#OK to Post
OR:
#Please Do Not Post
Arbi(S0,K,r,Sd,Su,T,P,x1,x2)
and write a new procedure
Arbi1(S0,K,r,Sd,Su,T,P,x2)
that does not input x1, but inside the program x1 is set to P-x2*S0. In other words, the seller of the option uses the full amount he or she got from the buyer of the option to build his portfolio.
Goodx2(S0,K,r,Sd,Su,T,P,h)
that inputs the same as above (except for x2) and a mesh-size h, and by trying out all x2 from 0 to 1, in increments of h, finds the set of x2's (that happens to be an interval), that allows arbitrage. Of course, it could be the empty set (interval). What are
EmpiOptionPrice(S0,K,r,Sd,Su,T,h1,h2)
that inputs everything as above, except for P, and a mesh sizes h1, h2, tries out all the P's, starting at 0, with increment h2, examines whether Goodx2(S0,K,r,Sd,Su,T,P,h1) is empty, and stops as soon as it is not empty, and returns the largest P for which it is empty.
What are
Contains procedures:
#Name
#Date, Assignment hw12.txt
#OK to Post
OR:
#Please Do Not Post
NOTE: IT IS STRONGLY RECOMMENDED THAT YOU DO THE HOMEWORK BY NEXT CLASS (March 5)
OP2scr(S1,S2,D11,D12,D13,D21,D22,D23,r,T,C1,C2,C3)
For the FAIR price (no arbitrage!) of an option that pays C1 dollars if the economy is bad, C2 dollars if the economy is OK, and C3 dollars if the economy is great, assuming inerest rate r.
Hint: you are allowed to construct portfolios consisting of so much Goverment bonds, so much from stock S1 and so much for stock S2 (so that you have three degrees of `freedom', in order to meet the three conditions that if the economy is bad you would give the buyer C1, if the economy is OK, you would give C2, and if the economy is good, you would give C3. You solve a system of three linear equations and three unknowns, getting the portfolio, and the fair price is today's price of that portfolio.
OP2nice(S1,S2,D11,D12,D13,D21,D22,D23,r,T,C1,C2,C3)
that outputs the same thing as OP2scr(S1,S2,D11,D12,D13,D21,D22,D23,r,T,C1,C2,C3)
Make sure that you do get the same thing!
Contains procedures:
OptionPrice(S0,D1,C): Inputs a list S0, or length N, say, where S0[i] is today's price of stock number i, D1 and C as above. Outputs the fair (no arbitrage) price of an option described by the list C, of length n, where C[j] is the seller's commitment to pay the buyer if the economy is in state j, and the list of lists D1 of length N, is such that D1[i][j] is the future price of stock number i IF the economy would wind-up in state number j
#Name
#Date, Assignment hw13.txt
#OK to Post
OR:
#Please Do Not Post
At time 0,
FindStateVector(S0,D1):
that inputs a list,S0, of current stock prices (as above), a list of lists of possible stock prices, D1, (with n=N) and outputs a state-price vector if it exists, and otherwise returns FAIL. What is the state-price-vector in the above data?
ArrowDebreu(D1,i)
That inputs a list of lists D1 (with n=N) as above, and an integer i between 1 and n, and outputs the ith Arrow-Debreu security defined on p. 11.
What are the Arrow-Debreu four securities for scenarios A,B,C,D, with the above data?
OptionPrice(S0,D1,C)
with S0=[1,s0], D1=[[exp(r*T),exp(r*T)],[su,sd]], and the needed C (you figure it out!), to prove (rigorously!) Exercise 8 (with t=0), on p. 19 of the book, where K is the strike price.
Contains, in addition to the procedures of C13.txt (use Help13(); for a list of them) procedures:
#Name
#Date, Assignment hw14.txt
#OK to Post
OR:
#Please Do Not Post
PrC(d,u,k) and PrP(d,u,k) ,
for the prices of a European Call option, and Put option, respectively, with strike price k (where in applications k is between d and u)
[Recall that for a Call option the seller has to give the buyer u-k if the stock went up, and nothing if the stock went down,
and for a Put option the seller has to give the buyer nothing if the stock went up, and k-d if the stock went down.]
Write four procedures
ECB(d,u,k,q), ECS(d,u,k,q), EPB(d,u,k,q), EPS(d,u,k,q),
For the expected gain
respectively, if they believe that the probability of the stock goind up is q (and hence, that the probability that the stock goes down is 1-q.)
PriceEO(CER,FERD,FERU,u,r,T,Pou,Dol)
That inputs
Check your procedure with the specific numbers in the example, that only has r,u, and T, symbolic.
Hint: the pdf is defined from 0 to infinity. The i-th moment (NOT about the mean) of a continuous random variable given by a pdf f(x) is the integral of xi f(x) over the domain. You also have to tell Maple that sigma is positive, so please do
assume(sigma,positive) ;
Another Hint: if m1, m2, m3, m4, etc., are the usual moments, then
Contains procedures:
#Name
#Date, Assignment hw15.txt
#OK to Post
OR:
#Please Do Not Post
Happy(w)
that inputs a sequence in {-1,1} and outputs the number of places where the partial sums are positive, or zero, but the previous partial sum is positive.
P(n,t)
that inputs a positive integer n, and outputs the sum of tHappy(w) over all sequences with n 1's and n -1's.
Can you conjecture an explicit expression for P(n,t)?
Q(n,t)
that inputs a positive integer n, and outputs the sum of tHappy(w) over all sequences in {1,-1} of length n (i.e. AllWalks(n,0) union AllWalks(n-1,1) ... union AllWalks(0,n))
Can you conjecture an explicit expression for the coefficients of Q(n,t)?
Contains procedures:
#Name
#Date, Assignment hw16.txt
#OK to Post
OR:
#Please Do Not Post
Use "backwards induction", and procedure FL (many times!), to write a procedure
PO(H,r,dt,C)
that inputs all possible stock-price histories, where H is a list of size, k+1, say, where H[i+1][j+1] is the price of the stock at time t=i if the history of downs and ups (starting at time 0) is coded by the binary representation of the integer j (j is between 0 meaning it was always down before time t=i, and 2^i-1, meaning it was always up), r is interest rate and dt is the duration of one time-unit. It also inputs the list of length 2^k-1 , where C[j+1] is the value of the claim if the history up to time t=k was coded by the integer j [j between 0 and 2^k-1, of course].
It outputs a list FC, analogous to H where FC[k+1] is C, and
FC[i+1][j+1] is the value of the option at time t=i if the history, up to that time
is coded by j. In particular, FC[1] should be a list of length 1, whose
value is the fair price of the option.
Test your procedure by a few random H's generated by procedure RandH.
Also with the 3-generation trees given in Fig. 2.3 and 2.4 (call and put options
resepectively) from the book.
In addition to the previous procedures from C16.txt, it contains procedures:
#Name
#Date, Assignment hw17.txt
#OK to Post
OR:
#Please Do Not Post
DrawBinTg(T,r,dt,C,eps)
that in addition to the Binomial tree, T, also inputs an interest rate r, and time unit dt (but when you try them out make r=0) a list the same length of nops(T[nops(T)]) that indicate the values of the claims corresponding to the stock price, and eps as before (eps=.2 looks good) and outputs a drawing of the tree, that in addition to indicating the stock prices at every node in black, also indicates the current value of the claim in red, and phi, (what portion of the stock to hold at that node), in blue.
Hint: use the option color=red etc. in the Maple plot command.
Contains procedures:
#Name
#Date, Assignment hw18.txt
#OK to Post
OR:
#Please Do Not Post
RandMar(k,M,U,K)
that we started in class, and I finished later, and understand the format of its output. Use it a few times with various k, M,U,K, but don't make M too big.
RandDSP(k,M,U,K) ,
that 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.
[Added April 5, 2015: Aniket Shah pointed out hat I did not mention how to generate the random values of the Discrete Stochastic Process. Just use rand(1..1000) or whatever].
Write a procedure
IsMartingale(P)
that 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. Test it for random outputs of RandDSP and RandMar (and hopefully get true for the latter).
In addition to the procedures of C18.txt (do Help18(); for a list), it contains the new procedures:
[NuVertices,ListOfChildren,ListOfGenerations,q,X]
phi a "predictable" ( a function defined on internal vertices) and Z0, the value of the new martingale at 0 (the top) outputs the brand-new martingale with the same vertices and the same prob. distibution but different values, according to Eq. (2.6), p. 36 of Alison Etheridge's "A course in Financial Calculus".
#Name
#Date, Assignment hw19.txt
#OK to Post
OR:
#Please Do Not Post
Write a procedure
CheckNewMartA(k,M,U,K)
that first creates a random martingale, using, RandMar(k,M,U,K), then generates a random phi (a random list of numbers with the same length as the number of non-leaves of that martingale), then uses procedure NewMart with a random Z0, and then uses, procedure IsMartingale(P), that you hopefully did for hw18.txt, to see if the new martingale is indeed a martingale, as promised.
CheckNewMart(k,M,U,K,HowMany)
That runs the above procedure CheckNewMartA(k,M,U,K) HowMany times, and returns true if all of them are true, and false otherwise. Run
CheckNewMart(10,3,3,5,10) ;
and hopefully get "true".
[Added April 5, 2015: Melissa Romanus pointed out that HowMany=1000 took too long, so I replaced it by 10].
Findphi(M,X)
that inputs two martingales, M, and X (with the same underlying trees and prob. distribution (i.e. the same sample space and sequence of sigma fields), and outputs Z0 and phi. Test it by first applying NewMart with a random phi, and see if you can recapture the phi.
BSdiscreteCall(u,d,r,N,K), BSdiscretePut(u,d,r,N,K)
that computes, NUMERICALLY, the fair price of a stock whose initial price is 1 dollar, lasts N generations, and follows the BINOMIAL model, where the price of the stock, from one generation to the next, may either go up by a factor u (u>1), or down with a factor d (d<1), with interest rate r, and strike price K.
What are
BSdiscreteCall(2,1/2,0,30,1.1), BSdiscretePut(2,1/2,0,30,1.1) ?
deriving the famous Black-Scholes formula.
Contains procedures:
#Name
#Date, Assignment hw20.txt
#OK to Post
OR:
#Please Do Not Post
Replace T by T-t in BSsym(S,sig,r,T,K), and call it V(S,t). use Maple to prove that V(S,t) satisfies the famous Black-Scholes PDE as given in wikipedia.
Added April 9, 2015: See the students' nice Solutions to the 20th homework assignment
Starting Brownian Motion Done right.
#Name
#Date, Assignment hw21.txt
#OK to Post
OR:
#Please Do Not Post
Added April 13, 2015: See the students' nice Solutions to the 21st homework assignment
#Name
#Date, Assignment hw22.txt
#OK to Post
OR:
#Please Do Not Post
[Version of April 14, 2015, 12:12pm, thanks to Anthony Zaleski]
and write a procdure
MyTaylor(f,x,a,N)
That inputs an expression f in the variable x, a number (or symbol) a, and a positive integer N, and outputs the polynomial of degree N consisting of the first N terms of the Taylor expansion of f at x=a.
Check the output of your program with the Maple taylor command for f=exp(x^2+x),a=0, and f(x)=sec(x), a=Pi.
Taylor2Start(f,x,y,h,k,N)
that inputs an expression f in the variables x, and y, symbols h and k, and a positive integer N, and outputs all the terms of total degree <=N in the taylor expansion of f(x+h,y+k) around the point (x,y).
What is the output of
ExpWtr(t,r)
That uses Ito's formula, and a recursion (generalizing example 4.3.2 on p. 87) to compute the Expectation of the Stochastic process
(Wt)r
for any specific (numeric) NON-NEGATIVE integer r, yielding an expression in t. Can you conjecture, and then PROVE a closed form (symbolic) expression (in r and t) for that quantity?
Hint: The Expectation of Int(f(W_s)dW_s, s=0..t) is always zero, by the martingale property.
Ints(f,sB,,B0,t0,N)
that inputs an expression, f in s and B, B0, a specific discete sample Brownian motion (gotten from Bt(t0,N)), a time t0, and a parameter N (a perfect square) and outputs a discrete approximation to the usual (classical) integral of f(s,B_0) with respect to s from 0 to t0.
IntW(f,s,B,B0,t0,N)
that inputs an expression, f in s and B, and also B0, a specific discete sample Brownian motion (gotten from Bt(t0,N), a time t0, and a parameter N (a perfect square) and outputs a discrete approximation to the STOCHASTIC integral of f(s,B_0) with respect to dW_s, from 0 to t0.
Added April 16, 2015: See the students' nice Solutions to the 22nd homework assignment
Contains procedures:
[Code by Anthony Zaleski]
[Code by Anthony Zaleski]
MakeSDE(f,t,Wt): inputs an expression f in t and Wt
and outputs an SDE , in Zt, whose solution is f(t,Wt). (Eq. 4.14) as a pair
[sig(t,Zt),mu(t,Zt)]
Meaning
dZt=sig(t,Zt) dWt + mu(t,Zt) dt
#Name
#Date, Assignment hw23.txt
#OK to Post
OR:
#Please Do Not Post
ExpftW(f,t,Wt,t0,N,M)
That inputs an expression f in t and Wt, a real number t0, a large integer N (a perfect square), and a fairly large integer M, that empirically estimates the expectation of the value of the Stochastic Process f(t,W_t) by taking M random discrete Brownian motions (using Bt(t0,N) M times), and averaging it.
Ints(f,s,B,B0,t0,N)
with M randomly-picked B0's from Bt(t0,N), write a procedure
ExpInts(f,s,B,t0,N,M)
that estimates the classical integral Int(f(s,B)ds)
IntW(f,s,B,B0,t0,N):
with M randomly-picked B0's from Bt(t0,N), write a procedure
ExpIntW(f,s,B,t0,N,M)
[Added April 18, 8:05pm: the previous version erroneusly had a B0 in the list of arguments (pointed out by Aniket Shat)]
that estimates the stochastic integral Int(f(s,W)dW)
CheckIto(f,s,Wt,t0,N,M)
that empirically (approximately, of course) verifies Ito's SDE for the Stochastic process f(t,W) .
f(t,Wt)
is given by the following integral formula
Int(f(t,x)*exp(-x^2/(2*t))/sqrt(2*Pi*t),x=-infinity..infinity)
Hint: prove it for monomials tmxn, using Maple's symbolic capacity, and the fact from hw22.txt about the expectation of powers of Wt^i . Then use the "linearity of expectation", i.e. using the fact that any function of x is a (usually "infinite") linear combination of monomials.
Read and understand the essay about Shepp's urn by W. Boyce.
Added April 20, 2015: See the students' nice Solutions to the 23rd homework assignment
Contains procedures:
#Name
#Date, Assignment hw24.txt
#OK to Post
OR:
#Please Do Not Post
GainLS(m,p,S) ,
that only computes the final gain (not the duration)
and
SimGain(m,p,S,M) ,
that runs a Shepp Urn game with m minus balls and p plus balls, M times, and computes the average.
Neigbors(S)
that inputs a strategy, and outputs the set of strategies where only one component is changed, by either 1 or -1 (and the entries are never negative), so there are at most 2*nops(S) neighbors of a strategy S.
BestNeighbor(m,p,S,M)
that inputs positive integers m and p, a strategy, S, and a large integer M (about a 1000), and outputs, the neighboring strategy that gives the best score of SimGain(m,p,S',M) (among all the neighbors S' of S) if it is better than that of S, or returns S, if all the neighbors have a lower score.
BestStrat(m,p,M)
that starts with the "chicken" strategy S=[0$p] and by the above "hill-climbing" finds a locally best strategy for m minus balls and p plus balls.
What are
[Note: of course, since this is done by simulation, you may get different answers each time, but hopefully not too far from each other]
Added April 27, 2015: See the students' nice Solutions to the 24th homework assignment
Contains procedures:
#Name
#Date, Assignment hw25.txt
#OK to Post
OR:
#Please Do Not Post
Vg(m,p,MaxLoss)
that computes the expected gain where the player refuses to ever have m-p more than MaxLoss (in other words, he or she do not want to take a chance, however small, of losing more than MaxLoss dollars).
What are the values of the sequences
evalf([seq(Vg(i,i,MaxLoss),i=1..100)]), for MaxLoss from 1 to 5?
How do they compare with evalf([seq(V(i,i),i=1..100)])?
Vxg(m,p,x,MaxLoss)
that finds the probability generating funcion with the above restiction
GuessPol1(L,n,d)
that inputs a list L, a symbol n, and a positive integer d (<=nops(L)-3), and outputs a polynomial of degree d in n, let's call it POL, such that L[i]=POL(i) (for i=1,2, ..., nops(L)). If it fails, it returns FAIL.
[Hint added, April 25, 2015, 7:10pm: create a generic polynomial, with undetermined coefficients, POL=a0+a1*n+...+ad*n^d, and create a set of unknowns, call them var, {a0,a1,...,ad}. Then create a set of nops(L) equations, call them, eq, {subs(n=1,POL)=L[1], ..., subs(n=nops(L),POL)=L[nops(L)]}. Then use Maple's solve command, to solve this system of linear equation. If the output of solve(eq,var) is NULL then return FAIL, otherwise, use the subs command to incorporate the output into the formerly generic POL, to get a concrete answer, if it exists].
GuessPol(L,n)
that inputs a list of numbers L, and outputs a polynomial of degree <=nops(L)-3, if it exits, or returns FAIL.
BdP(p,m0)
that inputs a SYMBOL p, a positive integer m0, and outputs the smallest integer, p0, such that Bd(p0,m0) is non-zero, and conjectures an explicit polynomial expressions for Bd(p,m0) (valid for p>=p0). Find them for m0 from 0 to 20.
ProveBdP(p,M0)
that first uses the above BdP to GUESS polynomial equations for Bd(p,m) and m0<=M0, and then goes on to prove these guesses rigorously (by induction), by using the recurrence (10) in William M. Boyce's article On a simple optimal stopping problem.
Added April 27, 2015: See the students' nice Solutions to the 25th homework assignment
Contains procedures:
#Name
#Date, Assignment hw26.txt
#OK to Post
OR:
#Please Do Not Post
[Corrected, Wed. April 29, 7:50am, thanks for Aniket Shah; Corrected again, April 30, 2015, thanks to John Kim]
Adapt procedure direct-pi, given in Werner Krauth's wonderful book (Algorihm 1.1) to write a Maple procedure
EstimateArea(F,x,y,N)
that inputs a polynomial F in x,y with positive coefficients, and a large positive integer N, and estimates the area of the region
{(x,y); F(x,y)2<1, x>=0, y>=0}
by picking N random points in a rectangle that contains the above region , and counts those that are inside it.
[Hint: first determine a rectangle containing the region, as small as possible (it is OK to do it by experimentation]
What did you get for
Adapt procedure markov-pi, given in Werner Krauth's wonderful book (Algorihm 1.2) to write a Maple procedure
EstimateAreaMarkov(F,x,y,N,delta)
that inputs a polynomial F in x,y with positive coefficients, and a large positive integer N, and a delta less than 1, and estimates the area of
{(x,y); F(x,y)2<1, x>=0, y>=0}
by picking N random points in a rectangle that contains the region F(x,y).
What did you get for
Added May 4, 2015: See the students' nice Solutions to the 26th homework assignment
Contains procedures:
#Name
#Date, Assignment hw27.txt
#OK to Post
OR:
#Please Do Not Post
OneStepIsing(M,z)
That inputs a matrix M (of size m by n, say) whose entries are {1,-1}, picks, uniformaly at random one location [i,j] (1<=i<=m, 1<=j<=n) and flips the M[i][j] entry (i.e. multiplies it by -1), getting a new matrix (that only differs from M in the [i,j] entry) definitely if Gibbs(M',z)/Gibbs(M,z) is >=1, and does the flipping with probability Gibbs(M',z)/Gibbs(M,z) if it is less 1 (using raR()). Otherwise it does nothing, and goes to the next step.
Ising(m,n,z,K)
that starts with an m by n {-1,1} matrix taken uniformly at random, a parameter z, and a very large positive integer K (say 30000), and only outputs the FINAL matrix.
DrawMat(M,h)
that enters a matrix (given a list of lists) with entries in {-1,1}, and "draws" it on a grid with resolution "pixel-size" h, where the pixel at location [i,j] is black if M[i][j]=1, and white if it is M[i][j]=-1
Ising(15,15,z,30000)
for z from 0 to 5 with increments of .05. Can you describe how it evolves?
Added May 4, 2015: See the students' nice Solutions to the 27th homework assignment
Contains procedures:
BEG(DI,K1): inputs a list DI of words and outputs the set of K1-prefixes
MagS1(DI,MAT): inputs a vocabulary list DI, and an already constructed matrix, MAT such that all the columns are STARTS of words in the vocabulary.
Renormalize(M)
that inputs a square matrix (given as a list of lists) with 3^k rows and 3^k columns (it should return FAIL if it is not of such a size) and outputs the 3^(k-1) by 3^(k-1) matrix obtained by splitting it into 3^(k-1) by 3^(k-1) blocks of 3 by 3 matrices, and replacing each such block by the SIGN of the sum of the entries (or if you wish, by "majority rule" if there are more 1's than -1's make the block 1, otherwise make it a -1)
Ising(81,81,z,5000)
for various z, and draw the shrinking matrices. See if you can get similar pictures to that Ken Wilson's article.
How many 3 by 3 "Magic Squares" are there (with E3). Can you find (and print) all of them?
[Added May 5, 2015, 6:55pm: if there are too many (i.e. if it takes more than 10 minutes, do instead the number of 3 by 3 magic squares only using the first 200 words of the vocabulary E3]
[Hint, to nicely print a matrix given as a list of lists, do "print(matrix(A));"]
[Added May 10, 2015: Matthew Kownacki, and Shalosh B. Ekhad, found the answer, see all 607904 3 by 3 magic squares
SMagS(DI,K1),
that outputs symmetric matices of letters, where the rows are identical to the corresponding columns.
FindSMS(DI,K1)
that finds just ONE symmetric K1 by K1 magic square using the vocabulary of DI. By using the file ENGLISH6, (where it is called E6), construct at least one such 6 by 6 symmetric magic square.
Added May 10, 2015: See the students' nice Solutions to the 28th homework assignment
Added May 18, 2014: Here are the students' evaluations.