**Overview**

MCMC algorithms used for simulating posterior distributions are indispensable tools in Bayesian analysis. A major consideration in MCMC simulations is that of convergence. Has the simulated Markov chain fully explored the target posterior distribution so far, or do we need longer simulations? A common approach in assessing MCMC convergence is based on running and analyzing the difference between multiple chains.

For a given Bayesian model, **bayesmh** is capable of producing multiple Markov chains with randomly dispersed initial values by using the **initrandom** option, available as of the update on 19 May 2016. In this post, I demonstrate the Gelman–Rubin diagnostic as a more formal test for convergence using multiple chains. For graphical diagnostics, see *Graphical diagnostics using multiple chains* in **[BAYES] bayesmh** for more details. To compute the Gelman–Rubin diagnostic, I use an unofficial command, **grubin**, which can be installed by typing the following in Stata: Read more…

**Initial thoughts**

Estimating causal relationships from data is one of the fundamental endeavors of researchers. Ideally, we could conduct a controlled experiment to estimate causal relations. However, conducting a controlled experiment may be infeasible. For example, education researchers cannot randomize education attainment and they must learn from observational data.

In the absence of experimental data, we construct models to capture the relevant features of the causal relationship we have an interest in, using observational data. Models are successful if the features we did not include can be ignored without affecting our ability to ascertain the causal relationship we are interested in. Sometimes, however, ignoring some features of reality results in models that yield relationships that cannot be interpreted causally. In a regression framework, depending on our discipline or our research question, we give a different name to this phenomenon: endogeneity, omitted confounders, omitted variable bias, simultaneity bias, selection bias, etc.

Below I show how we can understand many of these problems in a unified regression framework and use simulated data to illustrate how they affect estimation and inference. Read more…

\(

\newcommand{\xb}{{\bf x}}

\newcommand{\gb}{{\bf g}}

\newcommand{\Hb}{{\bf H}}

\newcommand{\Gb}{{\bf G}}

\newcommand{\Eb}{{\bf E}}

\newcommand{\betab}{\boldsymbol{\beta}}

\)I write ado-commands that estimate the parameters of an exponential conditional mean (ECM) model and a probit conditional mean (PCM) model by nonlinear least squares, using the methods that I discussed in the post Programming an estimation command in Stata: Nonlinear least-squares estimators. These commands will either share lots of code or repeat lots of code, because they are so similar. It is almost always better to share code than to repeat code. Shared code only needs to be changed in one place to add a feature or to fix a problem; repeated code must be changed everywhere. I introduce Mata libraries to share Mata functions across ado-commands, and I introduce wrapper commands to share ado-code.

This is the 27th post in the series **Programming an estimation command in Stata**. I recommend that you start at the beginning. See Programming an estimation command in Stata: A map to posted entries for a map to all the posts in this series.

**Ado-commands for ECM and PCM models**

I now convert the examples of Read more…

\(\newcommand{\xb}{{\bf x}}

\newcommand{\gb}{{\bf g}}

\newcommand{\Hb}{{\bf H}}

\newcommand{\Gb}{{\bf G}}

\newcommand{\Eb}{{\bf E}}

\newcommand{\betab}{\boldsymbol{\beta}}\)I want to write ado-commands to estimate the parameters of an exponential conditional mean (ECM) model and probit conditional mean (PCM) model by nonlinear least squares (NLS). Before I can write these commands, I need to show how to trick **optimize()** into performing the Gauss–Newton algorithm and apply this trick to these two problems.

This is the 26th post in the series **Programming an estimation command in Stata**. I recommend that you start at the beginning. See **Programming an estimation command in Stata: A map to posted entries** for a map to all the posts in this series.

**Gauss–Newton algorithm**

Gauss–Newton algorithms frequently perform better than Read more…

Autoregressive (AR) and moving-average (MA) models are combined to obtain ARMA models. The parameters of an ARMA model are typically estimated by maximizing a likelihood function assuming independently and identically distributed Gaussian errors. This is a rather strict assumption. If the underlying distribution of the error is nonnormal, does maximum likelihood estimation still work? The short answer is yes under certain regularity conditions and the estimator is known as the quasi-maximum likelihood estimator (QMLE) (White 1982).

In this post, I use Monte Carlo Simulations (MCS) to verify that the QMLE of a stationary and invertible ARMA model is consistent and asymptotically normal. See Yao and Brockwell (2006) for a formal proof. For an overview of performing MCS in Stata, refer to Monte Carlo simulations using Stata. Also see A simulation-based explanation of consistency and asymptotic normality for a discussion of performing such an exercise in Stata.

**Simulation**

Let’s begin by Read more…

**Overview**

In the frequentist approach to statistics, estimators are random variables because they are functions of random data. The finite-sample distributions of most of the estimators used in applied work are not known, because the estimators are complicated nonlinear functions of random data. These estimators have large-sample convergence properties that we use to approximate their behavior in finite samples.

Two key convergence properties are consistency and asymptotic normality. A consistent estimator gets arbitrarily close in probability to the true value. The distribution of an asymptotically normal estimator gets arbitrarily close to a normal distribution as the sample size increases. We use a recentered and rescaled version of this normal distribution to approximate the finite-sample distribution of our estimators.

I illustrate the meaning of consistency and asymptotic normality by Monte Carlo simulation (MCS). I use some of the Stata mechanics I discussed in Monte Carlo simulations using Stata.

**Consistent estimator**

A consistent estimator gets arbitrarily close in Read more…

\(\newcommand{\xb}{{\bf x}}

\newcommand{\betab}{\boldsymbol{\beta}}\)Before you use or distribute your estimation command, you should verify that it produces correct results and write a do-file that certifies that it does so. I discuss the processes of verifying and certifying an estimation command, and I present some techniques for writing a do-file that certifies **mypoisson5**, which I discussed in previous posts.

This is the twenty-fifth post in the series **Programming an estimation command in Stata**. I recommend that you start at the beginning. See Programming an estimation command in Stata: A map to posted entries for a map to all the posts in this series.

**Verification versus certification**

Verification is the process of establishing Read more…

This post was written jointly with Yulia Marchenko, Executive Director of Statistics, StataCorp.

As of update 03 Mar 2016, **bayesmh** provides a more convenient way of fitting distributions to the outcome variable. By design, **bayesmh** is a regression command, which models the mean of the outcome distribution as a function of predictors. There are cases when we do not have any predictors and want to model the outcome distribution directly. For example, we may want to fit a Poisson distribution or a binomial distribution to our outcome. This can now be done by specifying one of the four new distributions supported by **bayesmh** in the **likelihood()** option: **dexponential()**, **dbernoulli()**, **dbinomial()**, or **dpoisson()**. Previously, the suboption **noglmtransform** of **bayesmh**‘s option **likelihood()** was used to fit the exponential, binomial, and Poisson distributions to the outcome variable. This suboption continues to work but is now undocumented.

For examples, see *Beta-binomial model*, *Bayesian analysis of change-point problem*, and *Item response theory* under *Remarks and examples* in **[BAYES] bayesmh**.

We have also updated our earlier “Bayesian binary item response theory models using bayesmh” blog entry to use the new **dbernoulli()** specification when fitting 3PL, 4PL, and 5PL IRT models.

I make **predict** work after **mypoisson5** by writing an ado-command that computes the predictions and by having **mypoisson5** store the name of this new ado-command in **e(predict)**. The ado-command that computes predictions using the parameter estimates computed by ado-command **mytest** should be named **mytest_p**, by convention. In the next section, I discuss **mypoisson5_p**, which computes predictions after **mypoisson5**. In section Storing the name of the prediction command in e(predict), I show that storing the name **mypoisson5_p** in **e(predict)** requires only a one-line change to **mypoisson4.ado**, which I discussed in Programming an estimation command in Stata: Adding analytical derivatives to a poisson command using Mata.

This is the twenty-fourth post in the Read more…

**Overview**

I describe how to generate random numbers and discuss some features added in Stata 14. In particular, Stata 14 includes a new default random-number generator (RNG) called the Mersenne Twister (Matsumoto and Nishimura 1998), a new function that generates random integers, the ability to generate random numbers from an interval, and several new functions that generate random variates from nonuniform distributions.

**Random numbers from the uniform distribution**

In the example below, we use **runiform()** to create Read more…