> set.seed(1234)
> 
> r <- max(package)
> nT <- length(cases)
> 
> n <- rep(0,r)
> 
> n[1] <- length(cases[package==1])
> n[2] <- length(cases[package==2])
> n[3] <- length(cases[package==3])
> n[4] <- length(cases[package==4])
> 
> ybar <- rep(0,r)
> 
> ybar[1] <- mean(cases[1:n[1]])
> ybar[2] <- mean(cases[(n[1]+1):(n[1]+n[2])])
> ybar[3] <- mean(cases[(n[1]+n[2]+1):(n[1]+n[2]+n[3])])
> ybar[4] <- mean(cases[(n[1]+n[2]+n[3]+1):nT])
> ybar_all <- mean(cases)
> 
> SSTO <- sum((cases-ybar_all)^2)
> SSTR <- 0
> 
> for (i in 1:r) {
+ SSTR <- SSTR + n[i]*(ybar[i]-ybar_all)^2
+ }
> 
> SSE <- SSTO-SSTR
> 
> F_obs <- (SSTR/(r-1))/(SSE/(nT-r))
> 
> num_samp <- 10^5-1
> 
> F_obs_rand <- rep(0,num_samp)
> 
> 
> for (i1 in 1:num_samp) {
+ randgrp <- sample(nT,size=nT,replace=FALSE)
+ ybar[1] <- mean(cases[randgrp[1:n[1]]])
+ ybar[2] <- mean(cases[randgrp[(n[1]+1):(n[1]+n[2])]])
+ ybar[3] <- mean(cases[randgrp[(n[1]+n[2]+1):(n[1]+n[2]+n[3])]])
+ ybar[4] <- mean(cases[randgrp[(n[1]+n[2]+n[3]+1):nT]])
+ 
+ SSTR <- 0
+ 
+ for (i2 in 1:r) {
+ SSTR <- SSTR + n[i2]*(ybar[i2]-ybar_all)^2
+ }
+ 
+ SSE <- SSTO-SSTR
+ 
+ F_obs_rand[i1] <- (SSTR/(r-1))/(SSE/(nT-r))
+ }
> 
> (p_rand <- (sum(F_obs_rand >= F_obs)+1)/(num_samp+1))
[1] 5e-05
> 
> hist(F_obs_rand,xlab="F*",
+ main="Randomization Distribution for Cases Sold F-test")
> abline(v=F_obs)
>