

# Transition matrix  p(i,j) probability of landing on i from j
# 1-40 squares on the board, last number was not a double
# 41-80 last role was a double
# 81-120 last two rolls were doubles
# 121 just arrived in jail
# 122 one turn in jail
# 123 two turns in jail


# Initialise transition matrix
p=zeros(123);
#==============================================================================
# Dice probabilities
#==============================================================================
p1=zeros(12,1); # Excluding doubles
p2=zeros(12,1); # Doubles
for i=1:6
  for j=1:6
    if i!=j  
      p1(i+j)+=1/36;
    else
      p2(i+j)+=1/36;
    endif
  endfor
endfor
plot(1:12,p1,"r+",1:12,p2,"b+");
#==============================================================================
# Plain transitions: neglect special squares
#==============================================================================
# Outside jail
for i=1:40
  for j=1:12
    k=i+j-40*((i+j)>40);
    # Last throw wasn't a double
    p(   k,i)+=p1(j);
    p(40+k,i)+=p2(j);
    # Last throw was a double
    p(   k,40+i)+=p1(j);
    p(80+k,40+i)+=p2(j);
    # Last two throws were doubles
    p(   k,80+i)+=p1(j);
    p(121 ,80+i)+=p2(j); # Three doubles -> jail
  endfor
endfor
# Inside jail (square 11)
for j=1:12
  # Just arrived
  p(122  ,121)+=p1(j);
  p( 51+j,121)+=p2(j);
  # One turn in jail
  p(123  ,122)+=p1(j);
  p( 51+j,122)+=p2(j);
  # Last turn in jail have to come out
  p( 11+j,123)+=p1(j);
  p( 51+j,123)+=p2(j);
endfor
#==============================================================================
# Transition matrix for landing on "goto jail" square (no.31)
#==============================================================================
a=eye(123,123);
a( 31, 31)=0; a(121, 31)=1;
a( 71, 71)=0; a(121, 71)=1;
a(111,111)=0; a(121,111)=1;
p=a*p;
#==============================================================================
# Transition matrix for community chest, squares 3, 18 and 33 
#==============================================================================
a=eye(123,123);
for i=[3 18 33];
  for j=[0 40 80]
    a(  i+j,i+j)=14/16; # Remain in place
    a(  1+j,i+j)=1/16;  # Advance to go
    a(121  ,i+j)=1/16;  # Go (directly) to jail
    endfor
endfor
p=a*p;
#==============================================================================
# Transition matrix for chance, squares 8, 23, 37
#==============================================================================
a=eye(123,123);
for j=[0 40 80]
  for i=[8 23 37]
    a(   i+j,i+j)=6/16;  # Remain in place
    a(-3+i+j,i+j)+=1/16; # Go back three spaces
    a(   1+j,i+j)+=1/16; # Advance to Go
    a(   6+j,i+j)+=1/16; # Take a trip to Kings Cross
    a(  12+j,i+j)+=1/16; # Advance to Pall Mall
    a(  25+j,i+j)+=1/16; # Advance to Trafalger Square
    a(  40+j,i+j)+=1/16; # Advance to Mayfair
    a( 121  ,i+j)+=1/16; # Go (directly) to Jail
  endfor
  a(13+j, 8+j)+=1/16;
  a(29+j,23+j)+=1/16;
  a(13+j,37+j)+=1/16;  # Nearest utility 
  a(16+j, 8+j)+=2/16;
  a(26+j,23+j)+=2/16;
  a( 6+j,37+j)+=2/16;  # Nearest railway
endfor
p=a*p;
spy(p)
#==============================================================================
# Calculate probabilities
#==============================================================================
[v,lambda]=eig(p);
v1=v(:,1)/sum(v(:,1));
pp=v1(1:40)+v1(41:80)+v1(81:120);
pp(11)+=v1(121)+v1(122)+v1(123);
#==============================================================================
# Text plot
#==============================================================================
names=["Go"; "Old Kent Road"; "Community Chest"; "Whitechapel Road"; "Income Tax"; "Kings Cross Station"; "The Angel, Islington"; "Chance"; "Euston Road"; "Pentonville Road"; "Jail"; "Pall Mall"; "Electric Company"; "Whitehall"; "Northumberland Avenue"; "Marlybone Station"; "Bow Street"; "Community Chest"; "Marlborough Street"; "Vine Street"; "Free Parking"; "Strand"; "Chance"; "Fleet Street"; "Trafalger Square"; "Fenchurch St Station"; "Leicester Square"; "Coventry Street"; "Water Works"; "Piccadilly"; "Go to Jail"; "Regent Street"; "Oxford Street"; "Community Chest"; "Bond Street"; "Liverpool St Station"; "Chance"; "Park Lane"; "Super Tax"; "Mayfair"];
for i=1:length(pp)
  printf("%20s ",names(i,:))
  for j=1:floor(pp(i)*500)
    printf("*")
  endfor
  printf("\n")
endfor
#==============================================================================
# Top performing sets
#==============================================================================
set_names=["Railway"; "Brown"; "Light blue"; "Purple"; "Orange"; "Red"; "Yellow"; "Green"; "Bark Blue"]
value=zeros(9,1);
sigma=zeros(9,1);
value(1)=200*(pp(6)+pp(16)+pp(26)+pp(36));   
sigma(1)=200*sqrt(pp(6)+pp(16)+pp(26)+pp(36));       # Railways
value(2)=250*pp(2)+450*pp(4);
sigma(2)=sqrt(250^2*pp(2)+450^2*pp(4));              # Brown
value(3)=550*(pp(7)+pp(9))+600*pp(10);
sigma(3)=sqrt(550^2*(pp(7)+pp(9))+600^2*pp(10));     # Light Blue
value(4)=750*(pp(12)+pp(14))+900*pp(15);
sigma(4)=sqrt(750^2*(pp(12)+pp(14))+900^2*pp(15));   # Purple
value(5)=950*(pp(17)+pp(19))+1000*pp(20);
sigma(5)=sqrt(950^2*(pp(17)+pp(19))+1000^2*pp(20));  # Orange
value(6)=1050*(pp(22)+pp(24))+1100*pp(25);
sigma(6)=sqrt(1050^2*(pp(22)+pp(24))+1100^2*pp(25)); # Red
value(7)=1050*(pp(27)+pp(28))+1200*pp(30);
sigma(7)=sqrt(1050^2*(pp(27)+pp(28))+1200^2*pp(30)); # Yellow
value(8)=1275*(pp(32)+pp(33))+1400*pp(35);
sigma(8)=sqrt(1275^2*(pp(32)+pp(33))+1400^2*pp(35)); # Green
value(9)=1500*pp(38)+2000*pp(40);
sigma(9)=sqrt(1500^2*pp(38)+2000^2*pp(40));          # Dark blue






#==============================================================================
# Fair prices.
# What is a fair price for Old Kent Road, four players in the game,
# thirty minutes left on the clock.
# Player A owns Old Kent Road, Player B owns Whitechapel Road.
#==============================================================================
# Initial investment required by player b, 2 Hotels, 
initial_cost=500;
revenue_per_turn=250*pp(2)+450*pp(4)-230*(1/16)*(p(3)+p(18)+p(33))-200*(1/16)*(p(8)+p(23)+p(37));
meansq=250^2*pp(2)+450^2*pp(4)+230^2*(1/16)*(p(3)+p(18)+p(33))+200^2*(1/16)*(p(8)+p(23)+p(37));
number_of_turns=30/0.5; # 30 seconds per turn
profit=-initial_cost+revenue_per_turn*number_of_turns
# Market value of risk?
risk=sqrt(meansq*number_of_turns)
f_risk=0.1
fair_price=profit/2-f_risk*risk