Monte Carlo simulation is how we settle arguments in this course when the algebra gets hard to trust. You specify a complete, artificial data-generating process, use a random-number generator to draw a sample from that process, apply an econometric method, and repeat the exercise many times. Because you built the world, you know the true parameter and can see how the estimates behave around it. This does not replace a proof, but it gives us a numerical laboratory in which a theoretical claim has visible consequences.
Read the three uses on screen. First, we can confirm theory numerically. For example, theory says that OLS is unbiased when zero conditional mean and the other required assumptions hold. We can generate data that satisfy those assumptions and check whether repeated OLS estimates center on the coefficient we chose. Second, we can investigate a feature of real research data when no convenient theorem tells us what it will do. We reproduce that feature in a controlled simulation, change it while holding the rest of the process fixed, and watch the estimator respond. I use simulation this way in research as a diagnostic and design tool. Third, simulation helps you learn econometrics. A sampling distribution, bias, and variance become collections of actual numbers rather than only Greek letters. Today we begin with a result we already know, then deliberately violate an assumption and see the consequence.
It is a way to test econometric theories via simulation.
How is it used in econometrics?
Here is the question that motivates everything. In the regression model, zero conditional mean says that the average error is zero at every value of x. It rules out a systematic relationship between the observed regressor and the unobserved determinants collected in the error. We want to know what happens to OLS when that condition fails.
Could you simply take a real dataset and check? Open the answer and notice that both objects needed for a direct check are unobserved. You do not observe the true error, because the error includes all determinants of y that the model leaves out. An OLS residual is not a substitute for that true error. In a regression with an intercept, the OLS residuals are mechanically constructed to have zero sample covariance with x, even when the population error is related to x. You also do not observe the true coefficient beta. You observe beta-hat, one estimate from one sample, but without the truth you cannot label its gap as bias.
So real data hide both halves of the experiment: you cannot directly verify how the true error relates to x, and you cannot compare an estimate with the true parameter. Simulation restores both. We choose the coefficient and construct the error ourselves, then we can impose or violate the condition on purpose and measure exactly where repeated estimates center.
Suppose you are interested in checking what happens to OLS estimators if E[u|x]=0 (the error term and x are not correlated) is violated.
Question
Can you use the real data to do this?
Read the word in blue: you. You generate the data, and that single fact is what makes the whole method work. A data-generating process is a recipe specifying the parameter values, the distributions of the variables, the sample size, and any dependence among the variables. In real applications that recipe is unknown and must be inferred. Here, you write it down and make R follow it.
The first consequence on screen is that you know the true parameter. If you put a slope of one into the outcome equation, one remains the truth in every simulated sample. Each beta-hat will move because each sample is random, but you can compare the center of many beta-hats with one. That lets you distinguish ordinary sampling variation from systematic bias.
The second consequence is control. You can create two otherwise identical worlds and change only the feature you want to study, such as whether x and the error share a common component. Holding the coefficient, sample size, marginal distributions, and estimation command fixed prevents those other features from explaining the difference. This is a controlled experiment on an econometric method. In the first example we will make the assumptions hold and recover a familiar unbiasedness result. In the second, we will break one condition deliberately and watch the sampling distribution move.
You generate data (you have control over how data are generated)
Before we simulate anything, we need numbers that behave like random draws. The blue word “approximate” matters. A pseudorandom-number generator is a deterministic algorithm with an internal state. Starting from a state, it produces a long sequence whose frequencies, dependence, and other statistical properties are designed to resemble genuinely random draws. Physical sources of randomness exist, but pseudorandom generators are faster and, crucially for research, reproducible.
Run the first short cell. Say the call as “run-if five.” The argument five is the number of observations requested. Because the minimum and maximum arguments are omitted, R uses their defaults, zero and one. The result is five values from a Uniform zero-to-one distribution, so every interval of the same width inside zero to one has the same probability. The particular five values are not supposed to look evenly spaced. Randomness is a property of the process, not a promise about a tiny sample.
In the second cell, the assignment arrow stores ten thousand new Uniform draws in x, and “hist of x” groups them into bins. The horizontal axis records possible x values from zero to one. The vertical axis records the number of draws in each bin. Look for bars of roughly equal height across the range. They will not be exactly equal because ten thousand is still a finite random sample, but the approximately flat outline is the visible implication of a uniform density. The next tab explains why we can reproduce a particular sequence.
Pseudo random number generators (Pseudo RNG)
Algorithms for generating a sequence of numbers whose properties approximate the properties of sequences of random numbers
Examples
Draw from a uniform distribution:
Here is the practical consequence of “pseudo” rather than physically random. Once the generator’s algorithm and starting state are fixed, the entire future sequence is fixed. The first line says “set-seed two million three hundred eighty-seven thousand four hundred thirty-eight,” and that sets the starting state. The integer itself has no substantive meaning, and choosing a larger or more complicated seed does not make the simulation more random. It is simply a reproducible label for a position in the generator’s sequence.
Run the demonstration. The next call requests five Uniform zero-to-one draws, and with this seed the sequence begins approximately zero point zero four seven, zero point seven one two, zero point four zero seven, zero point two four two, and zero point three five seven. Run the whole cell again and you get the same five values, because the seed is reset before the draw. If you run only another random-draw command without resetting the seed, the generator advances and gives the next values in the sequence.
That answers the question at the bottom: setting a seed makes the entire analysis reproducible. A coauthor, a grader, or your future self can regenerate the same simulated data, estimates, tables, and figures. Put the seed once near the start of a simulation. Do not reset it inside every iteration, because that would reproduce the same sample over and over rather than generating independent replications. A seed does not repair a bad data-generating process, but it makes any problem traceable and any correct result auditable.
Numbers drawn using pseudo random number generators are not truly random
Demonstration
Question
What benefits does setting a seed have?
The Normal distribution is the workhorse for the examples that follow, and “R-norm” draws from it. Start with the left column. The mathematical notation says x follows Normal zero comma one, where zero is the mean and one is the variance. The call “R-norm ten thousand” asks for ten thousand observations. Because no other arguments are supplied, R uses mean zero and standard deviation one. The assignment stores those observations in x, and the histogram puts x values on the horizontal axis and bin counts on the vertical axis. You should see a bell-shaped distribution centered near zero, with most values relatively close to zero and fewer values in either tail.
Now look at the right column. The notation says Normal two comma four, so the population mean is two and the population variance is four. The R call gives the named mean argument the value two, which shifts the center from zero to two, and gives the named S-D argument the value two, which sets the standard deviation to two. This is the common trap: mathematical Normal notation on this slide reports variance as its second number, but R’s S-D argument requires the square root of that variance. Since the square root of four is two, the code is consistent with the notation. The right histogram should therefore be centered near two and spread about twice as widely in standard-deviation units as the left one. Small asymmetries and uneven bars are sampling noise, not a change in the distribution. In the next section, we organize R’s functions for the Normal and many other distributions.
x \sim N(0, 1)
x \sim N(2, 4)
This is a reference tab, and the naming pattern is more important than memorizing individual commands. The left column lists several distribution families built into R: Normal, Uniform, Beta, chi-square, F, logistic, and log-normal, plus many others. Each family has a short R name. Normal uses “norm,” Uniform uses “you-nif,” chi-square uses “chi-squared,” logistic uses “low-jiss,” and log-normal uses “ell-norm.”
The right column shows four prefixes that describe the operation you want. The “dee” prefix evaluates a density for a continuous distribution, or a probability mass for a discrete one. For example, dee-norm gives the height of a Normal density at a specified value. The “pee” prefix evaluates the cumulative distribution function, so pee-norm of q gives the probability that a Normal draw is less than or equal to q. The “cue” prefix reverses that mapping: give cue-norm a cumulative probability between zero and one, and it returns the corresponding cutoff value. Finally, the “R” prefix generates random draws, so R-norm of n returns n Normal observations.
You combine a prefix with a family name: dee-beta, pee-you-nif, cue-chi-squared, or R-low-jiss, for example. The distribution-specific arguments still matter, such as mean and standard deviation for a Normal, or minimum and maximum for a Uniform. Learn the four operations once and you can navigate a large part of R’s probability library. The next tabs use the standard Normal to make the distinction among density, cumulative probability, and quantile visible.
For each distribution, you have four different kinds of functions:
dnorm: density functionpnorm: distribution functionqnorm: quantile functionrnorm: random drawDee-norm of x gives the height of a Normal probability-density function at the value x. Because the mean and standard-deviation arguments are omitted here, both calls use the standard Normal, with mean zero and standard deviation one. The input is a location on the horizontal axis; the output is a density height on the vertical axis.
Read the figure from those axes. The red bell-shaped curve is the standard Normal density over values from minus three to three. The vertical line at zero marks the center, and the horizontal line marks density zero. A blue point at x equal to minus one is connected to the horizontal axis by a dashed guide. Its height is dee-norm of minus one, approximately zero point two four two. The other blue point is at x equal to two, where dee-norm of two is approximately zero point zero five four. Two is farther from the mean than minus one, so it lies lower in the tail. The Normal density is symmetric, which is why the height at minus one would equal the height at positive one, but not the height at two.
The crucial warning is that a density height is not itself a probability. For a continuous variable, the probability of drawing exactly one specified value is zero. Probabilities come from areas under the density over intervals. The units reinforce the distinction: a density can even exceed one for a sufficiently concentrated distribution, while a probability cannot. Flip to the pee-norm tab to see the area rather than the height.
dnorm(x) gives you the height of the density function at x.
dnorm(-1) and dnorm(2)
pnorm(x) gives you the probability that a single random draw is less than x.
Now the “pee” prefix, and this is the picture to keep in your head. Pee-norm of x is the Normal cumulative distribution function. With the default mean zero and standard deviation one, pee-norm of minus one asks for the probability that one standard Normal draw is less than or equal to minus one. For a continuous distribution, “less than” and “less than or equal to” have the same probability because the probability at the single endpoint is zero.
On the figure, the horizontal axis gives possible standard Normal values and the vertical axis gives density. The red curve is the same density from the preceding tab. The blue ribbon shades the area under that curve for x values at or below minus one. That area, not the curve’s height at minus one, is the requested probability. The annotation rounds the result to zero point one six; more precisely, R gives about zero point one five nine. So just under sixteen percent of standard Normal draws fall at least one standard deviation below the mean.
Compare the two functions directly. Dee-norm of minus one returned the height of the curve at minus one, about zero point two four two. Pee-norm of minus one accumulates every thin strip of area from the far left tail through minus one, giving about zero point one five nine. Density height and cumulative area answer different questions. The next tab moves the cutoff to the right side of the distribution but uses exactly the same cumulative operation.
This is the same cumulative-distribution function with the cutoff moved to two. The horizontal axis again shows standard Normal values, the vertical axis shows density, and the red curve is the standard Normal density. The blue ribbon now includes all area under the curve from the far left through x equal to two. Because two is two standard deviations above the mean, almost the whole bell is shaded.
Pee-norm of two is approximately zero point nine seven seven, which the figure rounds to zero point nine eight. Interpreted as a long-run proportion, about ninety-seven point seven percent of standard Normal draws fall below two. The unshaded right tail is the complement. Its probability is one minus pee-norm of two, approximately zero point zero two three, or two point three percent. R’s pee-norm function reports a left-tail cumulative probability by default, so subtraction from one is the standard way to form a right-tail probability.
Notice that dee-norm of two on the earlier tab was roughly zero point zero five four. That density height is not the two point three percent right-tail probability. The right-tail probability includes the entire area beyond two. This distinction becomes essential when we compute p-values: a test statistic locates a cutoff on the horizontal axis, while the p-value is tail area beyond that cutoff. The exercise tab now lets you use a nonstandard mean and standard deviation.
Your turn. The question asks for a cumulative probability under a Normal distribution, so use the “pee” prefix. In the work cell, enter pee-norm with a first argument of one, a named mean argument of one, and a named S-D argument of two. The first one is the cutoff whose left-tail probability you want. The named mean argument locates the center of the distribution at one, and the named S-D argument makes the standard deviation two, so the variance is four.
Before running the answer, standardize the cutoff. Subtract the mean from the value and divide by the standard deviation: one minus one, divided by two, equals zero. The requested cutoff is therefore zero standard deviations from the mean. A Normal distribution is symmetric around its mean, so exactly half of its area lies below the center and half lies above it. The answer is zero point five.
Open and run the answer cell to confirm. Read the code as “pee-norm of one, with mean one and S-D two.” Each argument is doing the job we just identified. The answer would remain zero point five for any positive standard deviation because changing the spread does not move the center. Do not generalize that irrelevance to other cutoffs. If the cutoff were not equal to the mean, changing the standard deviation would change how many standard deviations away it is and therefore change the probability.
What is the probability that a single random draw from a Normal distribution with mean = 1 and sd = 2 is less than 1?
Work here
sd isThe “cue” prefix reverses the question we have been asking. Pee-norm takes a cutoff on the measurement scale and returns the cumulative probability to its left. Cue-norm takes a cumulative probability and returns the cutoff that leaves that amount of area to its left. Read x in the definition as a probability, so it must be strictly between zero and one. Read q as the resulting value on the Normal variable’s scale.
For example, if cue-norm of x returns q, then the probability that a Normal draw falls below q is x. That is why the second bullet calls cue-norm the inverse of pee-norm: putting cue-norm of x inside pee-norm returns x, apart from tiny numerical rounding. The output is called the one-hundred-times-x percent quantile. A probability of zero point eight eight corresponds to the eighty-eighth percentile, the value with eighty-eight percent of the distribution below it.
As with the other Normal functions, the omitted defaults are mean zero and standard deviation one. So cue-norm of x by itself returns a standard Normal quantile. Supplying a different mean shifts the answer, and supplying a different standard deviation rescales its distance from the mean. Keep the direction straight: pee goes from value to probability; cue goes from probability to value. The next figure makes that inverse relationship concrete for zero point nine five.
qnorm(x), where 0 < x < 1, gives you the number q such that a single random draw is less than q with probability x.
pnorm(): pnorm(qnorm(x)) = xWe call the output of qnorm(x) the 100x\% quantile of the standard Normal distribution (because the default is mean = 0 and sd = 1 for qnorm()).
Here is a quantile you will meet repeatedly. The input to cue-norm is zero point nine five, so we want the point with ninety-five percent of the standard Normal distribution below it. R returns approximately one point six four four nine, usually rounded to one point six four or one point six four five.
Read the figure from the axes. The horizontal axis is the standard Normal value, and the vertical axis is density. The red line traces the density. The green ribbon shades the area for values below about one point six four, and that shaded area is ninety-five percent of the total area under the curve. The black point on the horizontal axis marks the cutoff, while the two annotations name the command and the rounded result. The remaining unshaded right tail has probability five percent.
This is the one-sided five-percent critical value. If a one-sided upper-tail test uses a standard Normal reference distribution, only five percent of null draws exceed about one point six four five. For a two-sided five-percent test, the five percent is split into two tails of two point five percent each. The positive cutoff is then cue-norm of zero point nine seven five, about one point nine six, with a matching negative cutoff. Those values are not magic constants. They are Normal quantiles chosen to leave specified tail areas, and the same probability-to-cutoff operation will reappear in hypothesis tests and confidence intervals.
Your turn again. The phrase “eighty-eight percent quantile” tells you to use the “cue” function with zero point eight eight as the probability input. The distribution has mean zero and standard deviation nine, so read the direct call in the work cell as “cue-norm of zero point eight eight, with mean zero and S-D nine.” The probability argument selects the percentile, the mean argument sets the center, and the S-D argument sets the scale. Notice that nine is already a standard deviation here, not a variance.
There is also a useful shortcut. First find the standard Normal eighty-eighth percentile, cue-norm of zero point eight eight, which is approximately one point one seven five. A quantile from a Normal with mean mu and standard deviation sigma equals mu plus sigma times the corresponding standard Normal quantile. Spoken for this problem, that is zero plus nine times one point one seven five, giving approximately ten point five seven five. So eighty-eight percent of draws from this Normal distribution lie below about ten point five seven, and twelve percent lie above it.
Open the answer and run the direct call to confirm. If the mean were not zero, you would add that mean after scaling. This location-and-scale rule is why standardizing is so useful: subtract the mean and divide by the standard deviation to move to the standard Normal, use its probabilities or quantiles, and then reverse the transformation when you need the value in the original units.
What is the 88% quantile of Normal distribution with mean = 0 and sd = 9?
Work here
qnorm(0.88) is about 1.17 for the standard Normal, and scaling the standard deviation by 9 scales the quantile by 9 as wellHere is the recipe, and every Monte Carlo exercise in this course follows it. In step one, specify the data-generating process completely: choose the true parameter values, sample size, distributions of the variables, and any dependence among them. In step two, generate one artificial dataset from that process. In step three, apply the estimator or procedure you want to study, such as OLS or a sample mean, and save its result. In step four, return to the data-generation step and repeat with fresh random draws many times. In step five, summarize the saved estimates and compare them with the truth chosen in step one.
The question at the bottom asks why steps one through three are repeated. More precisely, the specification in step one stays fixed while data generation and estimation are repeated. An estimator is a random variable because a new random sample produces a new estimate. One sample gives one draw from its sampling distribution. That single value cannot reveal the distribution’s center, spread, skewness, or tail behavior, so it cannot tell us whether a miss is ordinary sampling noise or systematic bias.
Thousands of independent replications approximate the sampling distribution. Their average estimates the estimator’s expected value, their standard deviation estimates its standard error under the chosen process, and their histogram shows its shape. Increasing the number of replications, often called B, reduces Monte Carlo noise in those summaries at the usual square-root-of-B rate. It does not change the estimator itself or the sample size within a replication. Keep that distinction in mind when we build the first example one step at a time.
Question
Why do the steps 1-3 many many times?
Let’s start with something you already believe, so you can see the simulation machinery work before we point it at a surprising result. The sample mean is one over n times the sum of the observations. Calling it unbiased for the population mean means that if we repeatedly drew samples and recomputed that average, the center of all those sample means would equal the population expected value.
The equation on screen states this precisely. Take the expectation of one over n times the sum from i equals one through n of x sub i. By linearity of expectation, move the constant one over n outside and take the expectation of each term. Every x sub i is drawn from the same distribution, so every term has expected value E of x. Adding n copies and dividing by n leaves E of x. Independence is part of the sampling setup shown here and will help determine the variance of the mean, although linearity itself does not require independence.
Unbiasedness does not say that a particular sample mean equals the population mean. Almost every individual sample mean will differ. It says the positive and negative sampling errors balance in expectation across repetitions. We have already proved that algebraically; now we will generate Uniform data, whose mean we know, and check whether the simulated sampling distribution centers there. Confirming a known result is also a diagnostic: if this simulation fails, the problem is in our code or design before we ever reach Example 2.
Question
Is sample mean really an unbiased estimator of the expected value?
That is, is E[\frac{1}{n}\sum_{i=1}^n x_i] = E[x], where x_i is an independent random draw from the same distribution?
Here are steps one through three for a single replication. The comment headings beginning with hash marks are ignored by R; they label the parts of the simulation for us. Read the first executable line as “x gets run-if one hundred.” Run-if one hundred generates one hundred independent Uniform zero-to-one values, and the assignment arrow stores that sample in x. The sample size is therefore one hundred. For a Uniform distribution from zero to one, the population mean is the midpoint, zero point five, and the population variance is one twelfth. We know that truth because we chose the distribution.
The next line calls “mean of x,” which adds the one hundred observations and divides by one hundred. The result is assigned to the object named mean-x. The final line contains only that object name, so R prints the stored sample mean. Run the cell. The answer should be near zero point five, but it will almost certainly not equal zero point five. Run it again and run-if advances to a fresh sample, so the printed mean changes.
That movement is sampling variation, not a defect. Under this design, the standard deviation of a sample mean is the square root of one twelfth divided by one hundred, about zero point zero two nine. A gap of a few hundredths is therefore ordinary. One result cannot tell us whether the estimator is centered correctly. We need to keep the process fixed, draw many samples, and store each mean, which is the purpose of step four.
Step four is where the method earns its keep. The preceding cell produced one estimate from one random sample. Seeing that estimate above or below zero point five does not reveal whether the difference is bad luck or systematic bias. We therefore repeat the same data-generating and estimation steps many times, using fresh random numbers on each repetition, and save the resulting collection of estimates.
The blue word identifies the programming tool: a loop. A loop tells R to execute a block of commands once for each value in a sequence. In this simulation, each pass through the block will draw one hundred new Uniform values and compute one new sample mean. The parameter value, distribution, and within-sample size remain fixed. Only the random realization changes. That consistency is what makes the resulting collection an approximation to one well-defined sampling distribution.
We also need a place to store each result. Without storage, R could compute a thousand means and immediately discard each one, leaving nothing to compare with the truth. The upcoming code therefore combines three pieces: choose the number of replications, allocate a storage vector of that length, and use the loop counter to put each estimate into a different position. If you have not written a loop before, the next tab strips away the econometrics and shows only the repetition syntax. After that, we put the sample-mean code inside it.
This is a loop in isolation, doing nothing except printing its counter. The first assignment sets capital B equal to one thousand. B is our notation for the number of Monte Carlo replications. The expression on the loop line, spoken “one through B,” constructs the integer sequence one, two, three, and so on through one thousand.
Now read the loop header as “for each i in one through B.” On the first pass, R assigns one to i. It then executes every command between the opening and closing braces. Here the only command is “print i,” so R prints one. On the second pass, i becomes two and R prints two. This continues through i equal to nine hundred ninety-nine and then one thousand. After the body has run for the final value, the loop ends. The braces matter because they define the block to repeat; indentation helps humans see that structure but does not cause the repetition.
Run the cell and you will see one thousand printed scalar values scroll through the output. The words “print i” after the hash mark are a comment for the reader, and R ignores them. The counter i is useful for more than displaying progress. In the real simulation, it will identify both which replication is running and which position in a storage object should receive that replication’s estimate. Flip to the next tab and replace the print command with the complete data-generation and estimation steps.
R code
Verbally
For each of i in 1:B (1, 2, \dots, 1000), do print(i).
i takes the value of 1, and then print(1)i takes the value of 2, and then print(2)i takes the value of 999, and then print(999)i takes the value of 1000, and then print(1000)Now we have the real simulation loop. The cell option spoken “auto-run true” tells the browser to execute this cell automatically, so its stored results are available to the next tab. The first assignment sets B, the number of replications, to one thousand. This is not the sample size. Each replication will still contain one hundred observations.
The call “rep zero comma B” creates a vector containing B zeros, and the assignment gives it the descriptive name estimate-storage-mean. The zero is merely a placeholder, and B determines the vector’s length. Preallocating storage is both clearer and faster than repeatedly growing a vector inside a loop.
Now read the loop body. For each i from one through B, run-if one hundred generates a fresh sample of one hundred Uniform zero-to-one values and stores it in x. Mean of x computes that sample’s mean and stores it temporarily in mean-x. The last line says “estimate-storage-mean bracket i gets mean-x,” meaning the i-th position of the storage vector receives this replication’s mean. On iteration one, position one is replaced; on iteration two, position two is replaced; by iteration one thousand, every placeholder has been replaced by an estimate. The objects x and mean-x are overwritten on each pass, but no result is lost because it has already been filed in its indexed position.
Run the cell if it has not autorun. Nothing is printed because every executable line is an assignment or part of the loop. That silence is correct. The result is a vector of one thousand sample means in memory, ready for the numerical and graphical comparisons on the next tab.
Step five compares the estimates with the truth. Because x was drawn from a Uniform zero-to-one distribution, the true population mean is zero point five. The first cell applies the mean function to the entire vector estimate-storage-mean. This is not another sample mean within one dataset. It is the Monte Carlo average of one thousand separate sample-mean estimates, our numerical approximation to the estimator’s expected value. It should be close to zero point five. In a verified run it was about zero point four nine nine, and an error around a thousandth is typical because the Monte Carlo standard error of this average is roughly zero point zero two nine divided by the square root of one thousand.
The second cell keeps all one thousand estimates visible. Read its first line as “data-frame, x equals estimate-storage-mean.” That makes a data frame whose x column contains the estimates. Gee-gee-plot uses that data, and geom-histogram maps the estimates to the horizontal axis, divides their range into thirty bins, fills the bars gray, and draws white borders between bins. The vertical axis is the number of estimates in each bin. The title reminds you that these are one thousand sample means, not ten thousand raw Uniform observations.
The first vertical line has x-intercept zero point five, is blue, and marks the truth. The second has x-intercept equal to the Monte Carlo average, is red and dashed, and marks the center of the simulated estimates. They should sit almost on top of each other. The histogram should also look approximately bell-shaped because each estimate averages one hundred independent values. Its spread is about zero point zero two nine; in the verified run, individual means ranged roughly from zero point four zero to zero point five eight. Those individual misses do not contradict unbiasedness. Unbiasedness is a statement about the center across repeated samples, not equality in any one sample. Increasing B makes our picture and Monte Carlo average more stable. Increasing the within-sample size of one hundred would make the histogram itself narrower.
Compare your estimates with the true parameter
runif()Now we reach the question that motivated the lecture. The model in the callout says y equals beta zero plus beta one times x plus u. Beta zero is the population intercept, beta one is the population slope, x is the observed explanatory variable, and u collects the unobserved determinants of y. The hat in beta-one-hat matters: beta one is the fixed truth built into the data-generating process, while beta-one-hat is the random OLS estimate computed from one sample.
Zero conditional mean says that the average value of u is zero for every value of x. When it holds, x carries no systematic information about the omitted component of y, and under the other OLS conditions the slope estimator is centered on beta one. The callout asks what happens when that condition does not equal zero. We are not changing the true coefficient. We are changing the relationship between x and the error, then asking where repeated beta-one-hats center.
Real data do not let us run this direct experiment because the true u and true beta one are both unobserved. In simulation, we choose beta one, generate u, and control exactly how u depends on x. The next tab creates a shared random component that makes x and u move together. After estimating the model a thousand times, we can compare every estimated slope with the known value one and determine both the direction and size of the resulting bias.
Question
What happens to \hat{\beta}_1 if E[u|x]\ne 0 when estimating y=\beta_0+\beta_1 x + u?
Read the preparation first. B equals one thousand replications, N equals one hundred observations per replication, and “rep zero comma B” preallocates a vector with one position for each estimated slope. The simulation does not set a seed on this tab, so exact displayed values can change across runs, but the data-generating process and its theoretical implications do not.
Inside each loop iteration, R-norm of N creates N independent standard Normal draws called mu. The next R-norm of N creates another independent standard Normal component and adds mu to form x. A third independent standard Normal component is added to the same mu to form u. Sharing mu is the deliberate violation. When mu is high, both x and u tend to be high; when it is low, both tend to be low. Because all components are jointly Normal here, the conditional mean of u given x is one half times x, not zero.
The outcome line says y equals one plus x plus u. It therefore builds in a true intercept of one and a true slope on x of one. Read the next call as “data-frame, y equals y, x equals x.” It combines the two observed variables into a dataset with columns named y and x. The error u is intentionally not passed to the regression, just as it would not be observed in real data.
The fee-ols function from the fixest package fits OLS. Read the formula as “y tilde x,” meaning regress y on x with an intercept, and read the named data argument as telling the function where to find those variables. The fitted model is stored as reg. Its coefficient component is a named vector, and “reg dollar coefficient bracket x” selects the estimated coefficient on x rather than the intercept. Finally, “estimate-storage bracket i” stores that coefficient in the i-th position. Only the x-error relationship was broken; the coefficients, sample size, and estimator stay fixed, giving us the controlled comparison promised at the start.
Before looking at the output, predict it. A simulation is most informative when you work out what should happen and then use the result to test your understanding and code. Write x as an independent standard Normal component plus mu, and write u as a different independent standard Normal component plus the same mu. All three underlying components have variance one and are mutually independent.
Now expand the covariance of x and u. Cross-covariances between the independent components are zero. The only term shared by x and u is mu, so the surviving covariance is the variance of mu, which equals one. For the denominator, the variance of x is the variance of its private standard Normal component plus the variance of mu. Independence lets us add those variances, giving one plus one equals two.
The probability limit displayed at the bottom says beta-one-hat approaches the true beta one plus the covariance of x and u divided by the variance of x. The true slope was set to one. The covariance-to-variance ratio is one over two. Adding them gives one point five. Spoken economically, OLS attributes to x both the true one-unit effect and the positive movement in the omitted error that accompanies x. In this jointly Normal construction, the conditional mean of u given x is exactly one half x, so the regression function itself has slope one point five.
The prediction is therefore not merely “biased upward.” It is quantitative: the sampling distribution of the estimated slope should center near one point five rather than the true structural coefficient one. Write down one point five before moving on. The comparison tab checks whether the simulated estimates agree with that derivation.
Both x and u contain \mu, so they are correlated by construction:
And the bias of \hat{\beta}_1 is governed by exactly that ratio:
\hat{\beta}_1 \rightarrow \beta_1 + \frac{Cov(x, u)}{Var(x)} = 1 + \frac{1}{2} = 1.5
The first cell applies the mean function to the one thousand stored slope estimates. That Monte Carlo average should be close to one point five. In a verified run of this data-generating process it was about one point five zero three, consistent with the covariance calculation on the previous tab. Small departures from one point five are Monte Carlo noise.
The plot then shows the full collection. Read the first call as “data-frame, x equals estimate-storage”; it puts the estimates in one column. The histogram maps that column to the horizontal axis, uses thirty bins, fills them gray, and separates them with white borders. The horizontal axis is the estimate of beta one; the vertical axis is the number of replications in each bin. The title identifies the one thousand estimates. A solid blue vertical line at one marks the true structural slope. A solid red vertical line at one point five marks the predicted center. The estimates pile up around the red line, while the blue truth lies well outside the bulk.
Compare this picture with the one in lecture 01-3, where zero conditional mean held and the histogram centered on the truth. Here the estimator is precise around the wrong target. That distinction also explains the final bullet. Raising N, the number of observations in each simulated dataset, reduces sampling variance and makes the histogram narrower, but the covariance-to-variance term remains one half. The distribution therefore tightens around one point five, not one. More data can make you increasingly confident in the wrong number when the identifying condition fails. Raising B instead would only make our simulated picture of that distribution smoother; it would not make the estimator in any one sample more precise.
Unbiasedness of OLS estimator \rightarrow Seeing it), where E[u|x]=0 held and the histogram was centered on the truthThis optional example asks whether OLS estimates its own uncertainty correctly. The model at the top has outcome y, intercept beta zero, slope beta one, regressor x, and error u. The bullets specify a standard Normal regressor, a Normal error, and zero conditional mean, so x contains no information about the average error. Under homoskedastic errors with variance sigma squared, the conditional variance of beta-one-hat, holding the observed x values fixed, is sigma squared divided by S-S-T sub X.
The denominator S-S-T sub X means the sum, from i equals one through n, of x sub i minus x-bar, all squared. It measures the sample variation in x around its mean. More spread in x makes this denominator larger and the slope estimator more precise. Sigma squared is the population error variance. It is unknown in real data, so the second line replaces it with sigma-hat squared.
Look at that replacement carefully. For each observation, u-hat sub i is the OLS residual. Squaring and summing the residuals gives the residual sum of squares. Dividing by n minus two estimates sigma squared because this simple regression estimated two coefficients, an intercept and a slope, using two degrees of freedom. Multiplying that estimate by one over S-S-T sub X gives the estimated variance of beta-one-hat. A regression table’s standard error is the square root of this estimated variance.
The question “is it unbiased?” asks whether, with x held fixed and new errors repeatedly drawn, the average estimated variance equals the true conditional variance. We will again build a world where every ingredient is known, store the variance estimate from each regression, and compare its average and distribution with the formula.
Model
\begin{aligned} y = \beta_0 + \beta_1 x + u \\ \end{aligned}Variance of the OLS estimator
True Variance of \hat{\beta}_1: V(\hat{\beta}_1) = \frac{\sigma^2}{\sum_{i=1}^n (x_i-\bar{x})^2} = \frac{\sigma^2}{SST_X}
Its estimator: \widehat{V(\hat{\beta}_1)} =\frac{\hat{\sigma}^2}{SST_X} = \frac{\sum_{i=1}^n \hat{u}_i^2}{n-2} \times \frac{1}{SST_X}
Question
Does the estimator really work? (Is it unbiased?)
Start with reproducibility. The first line says “set-seed nine hundred three thousand four hundred seventy-eight.” It fixes the random-number sequence, so the browser calculation matches the hidden R calculation used to print S-S-T sub X on the next tab. B is ten thousand replications and N is one hundred observations per replication. The two calls spoken “rep zero comma B” preallocate one storage vector for slope estimates and another for estimated slope variances.
The x line sits before the loop, which is a deliberate design choice. R-norm of N draws one hundred standard Normal x values once, and those exact values are reused in all ten thousand regressions. Read the next expression as “sum of x minus mean of x, squared.” It subtracts x-bar from every x, squares the deviations, and adds them to calculate S-S-T sub X. Holding x fixed matches the conditional variance formula on the preceding tab. If we redrew x each time, the true conditional variance would change across iterations and one fixed comparison line would no longer answer the same question.
Inside the loop, R-norm of N generates a fresh standard Normal vector on every pass, and multiplying by two gives u standard deviation two and variance four. The next tab uses that implemented value of sigma squared. Read the outcome line as “y gets one plus x plus u”; it sets the true intercept and slope to one. The following data-frame call supplies observed columns y and x to fee-ols from the fixest package, whose formula, spoken “y tilde x,” fits OLS with an intercept.
The coefficient line selects the named x entry from the fitted coefficient vector and stores it in position i. The variance line says “vee-cov of reg, bracket x comma x,” which returns the estimated coefficient variance-covariance matrix and extracts its row-x, column-x diagonal element. That element is an estimated variance, not its square-root standard error. With the conventional independent-and-identically-distributed variance calculation used here, it equals residual sum of squares divided by n minus two, then divided by S-S-T sub X. Ten thousand repetitions make the Monte Carlo average and density much more stable than one thousand would, while the fixed sample size N remains one hundred.
Here is the numerical comparison. The fixed x draw gives S-S-T sub X of about one hundred twelve point zero seven. The simulation code made u equal to two times a standard Normal, so sigma squared is two squared, or four. Substituting those values into the conditional variance formula gives four divided by one hundred twelve point zero seven, approximately zero point zero three five seven.
The cell says “mean of V-beta-storage,” averaging the ten thousand diagonal variance estimates saved from vee-cov. With the displayed seed, the verified result is approximately zero point zero three five eight. That is extremely close to the true value zero point zero three five seven; the small gap is simulation noise. The comparison concerns the variance estimator. Taking square roots afterward produces standard errors, and a nonlinear square root need not itself be exactly unbiased even when the variance estimate is.
There is also an algebraic reason for the agreement. Conditional on the fixed x values and under the Normal homoskedastic model, the residual sum of squares divided by sigma squared has a chi-square distribution with n minus two degrees of freedom. Its expected value is n minus two. Dividing the residual sum of squares by n minus two therefore has expectation sigma squared. S-S-T sub X is fixed in this experiment, so dividing by it preserves unbiasedness. The simulation is confirming that result numerically.
This matters because conventional regression inference builds on the estimated coefficient variance. Its square root is the reported standard error; dividing a coefficient difference by that standard error gives a t-statistic; and confidence intervals use it to measure sampling uncertainty. Later we will change the variance formula when errors are heteroskedastic or dependent. Here, under the assumptions being simulated, the conventional estimator targets the correct conditional variance.
True Variance
V(\hat{\beta}) = 4/112.07 = 0.0357
Check
Your Estimates of Variance of \hat{\beta}_1?
One last look, because an average hides the sampling variation around it. Read the first call as “data-frame, x equals V-beta-storage.” It places the ten thousand estimated variances in a column named x. Gee-gee-plot passes that data to geom-density, which maps x to the horizontal axis and draws a smoothed density rather than histogram counts. The vertical axis is therefore density, scaled so the total area under the curve is one. Geom-v-line adds a reference line at four divided by S-S-T sub X, rounded to four decimal places, and theme-lecture applies the deck’s common figure styling.
The horizontal values are possible estimates of the variance of beta-one-hat. The true conditional variance is about zero point zero three five seven, so the reference line should pass near the center of the density. With n minus two equal to ninety-eight degrees of freedom, the estimated variance is proportional to a chi-square random variable. Its distribution is not exactly symmetric: it is bounded below by zero and its right tail extends farther than its left. The skewness is modest at ninety-eight degrees of freedom, but individual estimates can still land noticeably above or below the truth.
The average from the preceding tab lies very close to the true line, which is what unbiasedness promises. It does not promise that most estimates equal the truth, or that the single estimated variance in a real dataset is especially close. In real work you observe one draw from this distribution. Increasing B would only make this displayed density smoother. Increasing the actual sample size N would concentrate the estimator more tightly. Keep the same lesson from the sample-mean example: unbiased describes the center across repetitions, not accuracy in every sample.
This optional exercise isolates something the preceding variance formula has already suggested: how does variation in x affect the precision of the OLS slope? The screen shows two simple-regression data-generating processes with the same intercept beta zero, the same slope beta one, and the same error distribution. The only difference is the explanatory variable used in each equation.
X one follows Normal zero comma one, so it has mean zero, variance one, and standard deviation one. X two follows Normal zero comma nine, so it has mean zero, variance nine, and standard deviation three. The parenthetical emphasizes the standard-deviation comparison because R’s R-norm function will require one and three as its named S-D arguments, even though the mathematical notation displays variances one and nine. The error u is standard Normal in both equations. Zero conditional mean holds for both x one and x two, so neither slope estimator has a systematic relationship with its error. Both should therefore center on the same true slope.
Predict the difference before opening the solution. The conditional slope variance is sigma squared divided by the sum of squared x deviations. With the same N, x two’s variance is nine times x one’s, so its sum of squared deviations should be roughly nine times larger. Its slope variance should therefore be about one ninth as large, and its slope standard deviation, the square root of variance, about one third as large. More variation in x supplies more information about how y changes with x. The simulation will check that both estimators remain unbiased while their precision differs.
Using MC simulations, find out how the variation in x affects the OLS estimators
Model setup
\begin{align*} y = \beta_0 + \beta_1 x_1 + u \\ y = \beta_0 + \beta_1 x_2 + u \end{align*}The preparation sets B to one thousand replications and N to one hundred observations per equation. Read the storage call as “matrix zero comma B comma two.” It allocates an object with B rows and two columns: each row is one replication, column one will hold the slope from equation one, and column two the slope from equation two.
Within each loop pass, the first R-norm call has N as the number of draws and a named S-D argument of one, so it draws x one with standard deviation one. The second has N draws and S-D three, so it draws x two with standard deviation three and variance nine. A third R-norm of N draws the standard Normal error u. The same u vector enters both outcomes, making the two simulated experiments paired within an iteration, but each regressor is independently drawn and satisfies zero conditional mean. Read the first outcome line as “y one gets one plus x one plus u.” It and the corresponding y two line both set the true intercept and slope to one.
Data-table from the data.table package combines the four named columns into one table. The package name written before the double colon tells R exactly where to find that function. Each fee-ols call from fixest fits one OLS regression with an intercept: y one on x one, then y two on x two, using the shared table supplied through the named data argument. The coefficient component is named, so bracket x one or bracket x two selects the slope. Matrix indexing uses i for the row and one or two for the column, filing both estimates from each iteration side by side.
After the loop, say “estimate-storage bracket all rows comma column one.” The blank position before the comma selects every row, and the one selects column one; that vector is named beta-one-s. Column two similarly becomes beta-two-s. The two mean calls should both return values close to one, confirming that changing x variation did not create bias. The two S-D calls measure the spread of the simulated sampling distributions. The x-two estimates should have about one third the standard deviation of the x-one estimates. In a verified one-thousand-replication run, the ratio was about two point eight three rather than exactly three, an ordinary finite-simulation difference from the theoretical approximation. The final tab makes that precision difference visible.
Here is the same result as a picture. The first two lines create separate data tables. In each, the column x contains one equation’s stored coefficient estimates, and the column called type contains the label “Equation 1” or “Equation 2” for every row. R-bind stacks those tables into one long dataset so a single plot can distinguish the two groups. The figure-width and figure-height options set the drawing canvas to nine by six point five inches; they affect presentation size, not the simulation or estimates.
In gee-gee-plot, geom-density maps coefficient estimates to the horizontal axis and group labels to fill color. It draws one smoothed density per equation. The alpha argument equals zero point five, making each fill half-transparent so you can see the overlap. Scale-fill-discrete receives an empty name, which removes the legend title. X-lab names the horizontal axis “Coefficient Estimate,” the lecture theme supplies common styling, and the final theme call places the legend below the plot. The vertical axis is density, so each curve has total area one; compare centers and widths, not raw heights as counts.
Both densities should center near the true slope of one. That is the unbiasedness result and is not what separates them. Equation one, whose x has variance one, produces the wider density. Equation two, whose x has variance nine, produces a dramatically narrower density, with roughly one third the standard deviation and one ninth the variance. The sample size, error variance, coefficient, and identifying condition are the same, so the difference comes from the denominator S-S-T sub X. More variation in an independent variable is information about its slope. If x barely moves, OLS has little leverage for separating changes in y associated with x from noise. A larger sample can also improve precision, but here we held N fixed to isolate the role of x variation.