R Asymptotics - Revision history

Rb: /* Analyse the distribution for estimated coefficients */

2016-02-02T22:57:03Z

‎Analyse the distribution for estimated coefficients

Rb: /* Analyse the distribution for estimated coefficients */

2016-02-02T22:56:20Z

‎Analyse the distribution for estimated coefficients

Rb: /* Draw the samples */

2016-02-02T22:50:52Z

‎Draw the samples

Rb at 22:50, 2 February 2016

2016-02-02T22:50:09Z

Show changes

Rb at 16:38, 6 February 2015

2015-02-06T16:38:44Z

Rb at 16:37, 6 February 2015

2015-02-06T16:37:55Z

Rb at 08:34, 6 February 2015

2015-02-06T08:34:57Z

Rb at 08:34, 6 February 2015

2015-02-06T08:34:10Z

Rb at 08:32, 6 February 2015

2015-02-06T08:32:31Z

Rb: Created page with "In this code we draw samples of different size from a population of data. We then apply OLS to each of these samples and obtain OLS parameter estimate distributions. We can sh..."

2015-02-06T08:31:22Z

Created page with "In this code we draw samples of different size from a population of data. We then apply OLS to each of these samples and obtain OLS parameter estimate distributions. We can sh..."

New page

In this code we draw samples of different size from a population of data. We then apply OLS to each of these samples and obtain OLS parameter estimate distributions. We can show how the parameter estimate variance reduces with increased sample sizes.

# Code to demonstrate unbiasedness of OLS estimator
# also the bahviour as sample size increases

Npop <- 100000 # Population Size
Nsamp1 <- 1000 # Sample Size
Nsamp2 <- 10000 # Sample Size
Nsamp3 <- 100000 # Sample Size

Q <- 600 # number of samples taken

u <- rnorm(Npop) # true error terms
beta <- matrix(c(0.5, 1.5),2,1) # true parameters
Xpop <- matrix(1,Npop,2) # initialise population X
Xpop[,2] <- floor(runif(Npop)*12+8) # exp var uniform r.v. in [8,19]
Ypop <- Xpop %*% beta + u # %*% is matrix multiplication

save_beta1 <- matrix(0,Q,3)

for (i in 1:Q ) {
sel1 <- ceiling(runif(Nsamp1)*Npop) # select Nsamp1 observations
Ysel <- Ypop[sel1,]
Xsel <- Xpop[sel1,]
reg_sel <- lm(Ysel~Xsel[,2]) # estimate regression in sample
save_beta1[i,1] <- reg_sel$coefficients[2] # save estimated beta_1

sel2 <- ceiling(runif(Nsamp2)*Npop) # select Nsamp2 observations
Ysel <- Ypop[sel2,]
Xsel <- Xpop[sel2,]
reg_sel <- lm(Ysel~Xsel[,2]) # estimate regression in sample
save_beta1[i,2] <- reg_sel$coefficients[2] # save estimated beta_1

sel3 <- ceiling(runif(Nsamp3)*Npop) # select Nsamp3 observations
Ysel <- Ypop[sel3,]
Xsel <- Xpop[sel3,]
reg_sel <- lm(Ysel~Xsel[,2]) # estimate regression in sample
save_beta1[i,3] <- reg_sel$coefficients[2] # save estimated beta_1
}

# define bins
bins = seq(min(save_beta1[,1])-0.05, max(save_beta1[,1])+0.05, by = 0.005) # match

# Kernel Density Estimates (smooth histograms)
d1 <- density(save_beta1[,1])
d2 <- density(save_beta1[,2])
d3 <- density(save_beta1[,3])

# Plots
colors <- rainbow(3)
plot(range(c(1.35,1.65)), range(c(0,max(d3$y)+0.05)), type="n", xlab="beta", ylab="Density")
# axis(1, at = c(1,2,3,4,5), labels=c("Fail","3", "2.2", "2.1", "1"))
lines(d1$x,d1$y, col = colors[1], lwd=3)
lines(d2$x,d2$y, col = colors[2], lwd=3)
lines(d3$x,d3$y, col = colors[3], lwd=3)
title(main = "Distribution of estimated Parameters")
legend("topleft", inset=.0, title="Sample Size",c("Nsamp1","Nsamp2","Nsamp3"), fill=rainbow(3), bty = "n")

@@ Line 75: / Line 75: @@
     legend(&quot;topleft&quot;, inset=.0, title=&quot;Sample Size&quot;,
            c(&quot;Nsamp1&quot;,&quot;Nsamp2&quot;,&quot;Nsamp3&quot;), fill=rainbow(3), bty = &quot;n&quot;)</pre>
-[[File:OLSasympdistr_kernel.jpeg|frame|none]]
+[[File:OLSasympdist kernel.jpeg|frame|none]]

@@ Line 75: / Line 75: @@
     legend(&quot;topleft&quot;, inset=.0, title=&quot;Sample Size&quot;,
            c(&quot;Nsamp1&quot;,&quot;Nsamp2&quot;,&quot;Nsamp3&quot;), fill=rainbow(3), bty = &quot;n&quot;)</pre>
-[[File:OLSasympdistr_files/figure-html/unnamed-chunk-5-1.png|frame|none]]
+[[File:OLSasympdistr_kernel.jpeg|frame|none]]

@@ Line 21: / Line 21: @@
 Now we have large vectors of <code>Ypop</code> and <code>Xpop</code> which contain the population data.
-= Draw the samples =
+= Prepare the sampling =
 We set up the number of samples to draw and the (three different) sample sizes
@@ Line 32: / Line 32: @@
 <pre class="r">save_beta1 &lt;- matrix(0,Q,3)</pre>
 = Now draw the samples =

@@ Line 51: / Line 51: @@
      plot(range(c(1.35,1.65)), range(c(0,max(d3<source enclose=none>$</source>y)+0.05)), type="n", xlab="beta", ylab="Density")
      # axis(1, at = c(1,2,3,4,5), labels=c("Fail","3", "2.2", "2.1", "1"))
-     lines(d1$x,d1$y, col = colors[1], lwd=3)
+     lines(d1<source enclose=none>$</source>x,d1<source enclose=none>$</source>y, col = colors[1], lwd=3)
-     lines(d2$x,d2$y, col = colors[2], lwd=3)
+     lines(d2<source enclose=none>$</source>x,d2<source enclose=none>$</source>y, col = colors[2], lwd=3)
-     lines(d3$x,d3$y, col = colors[3], lwd=3)
+     lines(d3<source enclose=none>$</source>x,d3<source enclose=none>$</source>y, col = colors[3], lwd=3)
      title(main = "Distribution of estimated Parameters")
      legend("topleft", inset=.0, title="Sample Size",c("Nsamp1","Nsamp2","Nsamp3"), fill=rainbow(3), bty = "n")

@@ Line 24: / Line 24: @@
        Xsel <- Xpop[sel1,]
        reg_sel <- lm(Ysel~Xsel[,2])   # estimate regression in sample
-       save_beta1[i,1] <- reg_sel\$coefficients[2]  # save estimated beta_1
+       save_beta1[i,1] <- reg_sel<source enclose=none>$</source>coefficients[2]  # save estimated beta_1
        sel2 <- ceiling(runif(Nsamp2)*Npop)    # select Nsamp2 observations

@@ Line 30: / Line 30: @@
        Xsel <- Xpop[sel2,]
        reg_sel <- lm(Ysel~Xsel[,2])   # estimate regression in sample
-       save_beta1[i,2] <- reg_sel$coefficients[2]  # save estimated beta_1
+       save_beta1[i,2] <- reg_sel<source enclose=none>$</source>coefficients[2]  # save estimated beta_1
        sel3 <- ceiling(runif(Nsamp3)*Npop)    # select Nsamp3 observations
@@ Line 36: / Line 36: @@
        Xsel <- Xpop[sel3,]
        reg_sel <- lm(Ysel~Xsel[,2])   # estimate regression in sample
-       save_beta1[i,3] <- reg_sel\$coefficients[2]  # save estimated beta_1
+       save_beta1[i,3] <- reg_sel<source enclose=none>$</source>coefficients[2]  # save estimated beta_1
+     }