################################
################################ St. Petersburg Paradox in 6 Lines of Code - Roland B. Stark, Integrative Statistics, June 2022
################################

# This is companion code to the piece at https://www.integrativestatistics.com/blog


# A.  Initial, most concrete and transparent, but most tedious version
 
N = 1000000 # Sets the number of trials
rand = runif (N,0,1) # A uniform random number, ranging from 0 to 1, to randomize results, with sample size matching our N above
pay = rep (2, N) # This initializes "pay" at the lowest possible value of $2, when the first flip is heads
pay [rand < 1/2] = 2^2 # Long-term, half the time the first flip will come up tails and the 2nd (or greater), heads, paying out at least $4
pay [rand < 1/4] = 2^3
pay [rand < 1/8] = 2^4
pay [rand < 1/16] = 2^5
pay [rand < 1/32] = 2^6
pay [rand < 1/64] = 2^7
pay [rand < 1/128] = 2^8
pay [rand < 1/256] = 2^9
pay [rand < 1/512] = 2^10
pay [rand < 1/2^10] = 2^11 # Starting here we simplify our task by using exponents instead of large numbers within the brackets
pay [rand < 1/2^11] = 2^12
pay [rand < 1/2^12] = 2^13
pay [rand < 1/2^13] = 2^14
pay [rand < 1/2^14] = 2^15
pay [rand < 1/2^15] = 2^16
pay [rand < 1/2^16] = 2^17
pay [rand < 1/2^17] = 2^18
pay [rand < 1/2^18] = 2^19
pay [rand < 1/2^19] = 2^20
pay [rand < 1/2^20] = 2^21
pay [rand < 1/2^21] = 2^22
pay [rand < 1/2^22] = 2^23
pay [rand < 1/2^23] = 2^24
pay [rand < 1/2^24] = 2^25
pay [rand < 1/2^25] = 2^26
pay [rand < 1/2^26] = 2^27
pay [rand < 1/2^27] = 2^28
pay [rand < 1/2^28] = 2^29
pay [rand < 1/2^29] = 2^30
pay [rand < 1/2^30] = 2^31
pay [rand < 1/2^31] = 2^32
pay [rand < 1/2^32] = 2^33
pay [rand < 1/2^33] = 2^34
pay [rand < 1/2^34] = 2^35 # 2^34 =~ 17 billion; an arbitrary stopping point

mean(pay); max(pay) 
# For me, with 8GB of RAM and a 2.6Ghz processor, this runs in 1 second for 1M trials
# The larger our N, the greater the chance of a very large payout among the trials, and the greater mean payout


################################

# B. Using a loop as a shortcut, we efficiently create the equivalent of the long set of commands in A. 

# For me, this loop runs in 1 second for 1M trials and about 20-40 seconds for 100M trials

N = 1000000
rand = runif (N,0,1)
pay = rep (2, N) 
for (i in 1:34) {  # This loops through to compute pay for each iteration of rand, from 2^1 through 2^34
pay [rand < (1/(2^(i)))] = 2^(i+1) } # e.g., when rand = 1/2^4, pay = 2^5
mean (pay); max (pay)


################################


# C.  Here's an alternative approach from George Pipis at https://predictivehacks.com/st-petersburg-paradox/

# For me, it takes much longer than the B. version to run:  1M trials required 31 minutes rather than 1 second


x = Sys.time() # Records the start time
payoffs=c()
 for (i in 1:1000000) {
   k=1
   while (sample(c("H","T"),1)!="H") {
    k=k+1
   }
   payoffs=c(payoffs, 2**k)
 }
 summary(payoffs)
Sys.time() - x # Gives the time elapsed since start time


################################


# D. SPSS code (created with SPSS v21)

*1st create large N of cases.

new file.
input program. 
loop #i=1 to 1000000. 
end case.
end loop. 
end file.
end input program.
compute temp = 1.
fre temp.
*This confirms that 1M cases have been created.

*Now for the main simulation.

do repeat x=v1 v2 v3 v4 v5 v6 v7 v8 v9 v10.
*One could use more variables.
compute x = rv.uniform(0,1).
recode x (lo thru .5 = 0) (.5 thru hi = 1).
end repeat.
exe.
*Now each variable, v1 through v10, consists of 0s and 1s, which can stand for tails and heads.

do if v1=1.
compute p=2.
else if v2=1.
compute p=4.
else if v3=1.
compute p=8.
else if v4=1.
compute p=16.
else if v5=1.
compute p=32.
else if v6=1.
compute p=64.
else if v7=1.
compute p=128.
else if v8=1.
compute p=256.
else if v9=1.
compute p=512.
else if v10=1.
compute p=1024.
end if.
exe.

var lab p 'Payout'.

graph/barchart p.
*This is a basic graph of payouts; a barchart is much more effective in this case than a histogram; see why:
graph/hist p.

*SPSS''s Chart Builder allows for a few more specs.
GGRAPH
  /GRAPHDATASET NAME="graphdataset" VARIABLES=p COUNT()[name="COUNT"] MISSING=LISTWISE 
    REPORTMISSING=NO
  /GRAPHSPEC SOURCE=INLINE.
BEGIN GPL
  SOURCE: s=userSource(id("graphdataset"))
  DATA: p=col(source(s), name("p"), unit.category())
  DATA: COUNT=col(source(s), name("COUNT"))
  GUIDE: axis(dim(1), label("Payout"))
  GUIDE: axis(dim(2), label("Percent"))
  SCALE: linear(dim(2), include(0))
  ELEMENT: interval(position(summary.percent(p*COUNT, base.all(acrossPanels()))), 
    shape.interior(shape.square))
END GPL.

*For me this whole set of SPSS code runs in 10 seconds.