# Coin Toss Simulation
outcomes <- c(0, 1) # Tails for 0 and Heads for 1
p <- 0.5 # true value of coming up Heads
nSims <- 10
Copyright By PowCoder代写 加微信 powcoder
meanFlips <- array(dim = nSims) for (i in 1:nSims) { flips <- sample(outcomes, size = T, replace = TRUE, prob = c(1-p, p)) meanFlips[i] <- mean(flips) # estimator for p # Each individual result is different from p # but on average, the estimator is unbiased # Demonstrating CLT outcomes <- c(1, 0) nSims <- 10000 Z <- array(dim = nSims) for (i in 1:nSims) { flips <- sample(outcomes, size = T, replace = TRUE, prob = c(p, 1-p)) Z[i] <- (mean(flips) - p) / sqrt(p * (1-p) / T) plot(density(Z), xlim = c(-3, 3)) curve(dnorm, -3, 3, add = TRUE, col = 'red') # Comparing t- and normal distribution # 2.5% & 97.5% quantile: qnorm(c(0.025, 0.975)) qt(c(0.025, 0.975), df = 5) qt(c(0.025, 0.975), df = 10) qt(c(0.025, 0.975), df = 15) qt(c(0.025, 0.975), df = 20) qt(c(0.025, 0.975), df = 25) qt(c(0.025, 0.975), df = 30) qt(c(0.025, 0.975), df = 50) qt(c(0.025, 0.975), df = 100) # Implementing the Jackknife procedure: library(quantmod) getSymbols('MSFT', from = '2007-01-01', to = '2021-12-31') MSFT <- Ad(MSFT) returns <- monthlyReturn(MSFT, type = 'log') T <- nrow(returns) # 99% CI using formulas for normal distribution approximation sigmaHat <- sd(returns) qnorm(c(0.005, 0.995)) c(sigmaHat - 2.58 * sigmaHat / sqrt(2*T), sigmaHat + 2.58 * sigmaHat / sqrt(2*T)) # Jackknife procedure muHat <- array(dim = T) for (i in 1:T) { muHat[i] <- mean(returns[-i]) muBar <- mean(muHat) # Jackknife mean (T-1) * (muBar - mean(returns)) # jackknife estimate of bias se.jack.muHat <- sqrt((T-1)/T * sum( (muHat - muBar)^2)) # jackknife estimate of se sd(returns) / sqrt(T) # actual standard error estimate from formula # histogram of the leave-one-out estimates: hist(muHat) # using bootstrap package install.packages('bootstrap') # use jackknife function from bootstrap package library(bootstrap) muHat.jack <- jackknife(coredata(returns), mean) muHat.jack$jack.se # The Bootstrap Procedure B <- 10000; muHat.boot <- array(dim = B) set.seed(100) for (i in 1:B) { boot.data <- sample(returns, T, replace = TRUE) # bootstrap sample muHat.boot[i] <- mean(boot.data) # bootstrap sample mean # bootstrap bias estimate: mean(muHat.boot) - mean(returns) # bootstrap estimate for standard error: sd(muHat.boot) # Original standard error sd(returns) / sqrt(T) # histogram of sample means hist(muHat.boot) qqnorm(muHat.boot) qqline(muHat.boot) # Bootstrap Percentiles quantile(muHat.boot, c(0.025, 0.975)) # Basic Method delta <- muHat.boot - mean(returns) c(mean(returns) - quantile(delta, 0.975), mean(returns) - quantile(delta, 0.025)) 程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com