## Example 1 

setwd("C:/local/classes/M483")
data <- read.csv('DJret.csv', header = T)
 names(data)
Y <- data$Ret.Pct.

Ytrunc <- pmax(pmin(Y,5),-5)

hist(Ytrunc, breaks = seq(-5,5,0.25), freq=F )

(Xbar <- mean(Y))
(s <- sd(Y))

 nd <- function(x){
   dnorm(x, Xbar, 0.85^2)
 }

 hist(Ytrunc, breaks = seq(-5,5,0.25), freq=F )
  curve(nd, add=T, col="red")

 td <- function(t){
   dt(t/a, df = df.t)/a
 }
 
 hist(Ytrunc, breaks = seq(-5,5,0.25), freq=F )
  df.t <- 5
  a <- 0.7
  curve(td, add=T, col="red", lwd=2)
 

qqnorm(Y)     ## normal QQ plot


## Example 2
## mean and median for n=3 from Uniform distribution

 M <- 10000    # simulation size
 n <- 3       # sample size

 X <- matrix(runif(n*M),M,n)
 
 means <- rep(0,M)
 medians <- rep(0,M)
 maxs <- matrix(0,M,1)

 for (i in 1:M){
   means[i] <- mean(X[i,])
   medians[i] <- median(X[i,])
   maxs[i] <- max(X[i,])
 }

#plot( density(means))
 hist(2*means, freq= FALSE)
  lines(density(2*medians), col="red", lwd=2)
  lines(density(maxs), col="blue", lwd=2)


M <- 5000    # simulation size
 n <- 30     # higher sample size

 X <- matrix(runif(n*M),M,n)
 
 means <- rep(0,M)
 maxs <- matrix(0,M,1)

 for (i in 1:M){
   means[i] <- mean(X[i,])
   maxs[i] <- max(X[i,])
 }

 hist(maxs, freq= FALSE, xlim=c(0.8,1.2))
  lines(density(2*means), col="blue", lwd=2)

  
## Example 3: heavy-tailed distribution 

 # Cauchy = t-dist with df=1 

  M <- 5000    # simulation size
  n <- 30     # higher sample size

 X <- matrix(rt(n*M, df = 1),M,n)

 means <- rep(0,M)
 medians <- matrix(0,M,1)

 for (i in 1:M){
   means[i] <- mean(X[i,])
   medians[i] <- median(X[i,])
 }

 hist(medians, freq= FALSE, xlim=c(-10,10))
  lines(density(means), col="blue", lwd=2)

  hist(means)