Simulating Data with MCMC with priors

I have data that are a sample of averages that each have a known variance/standard deviation. My goal is to have simulated data whose estimates vary in a similar way to the observed data. Through my own research, I was pointed to using Markov-Chain Monte Carlo.

My approach is to generate simulations of population forecasts that vary according to historical growth rates. The data that I am using are the natural log-transformed numbers of the forecasted population END and the variance END_VAR of the forecast. I want to generate samples/chains that are distributed with a mean of END and a variance of END_VAR.

The goal is to take the 50 estimates of the forecasted state populations in 2020 and compare them to each other for each simulation. Based on previous growth rates, some states grew faster during certain years such as the Great Recession (e.g. Texas) or lost population due to natural disasters such as Louisiana. I want to account for these historical growth rates with the standard deviation in the model. The end goal (after the MCMC step shown here) is to find the congressional apportionment for each state, but there is some amount of error in forecasting that population. Since a small difference in population for some states (e.g. California) can mean the difference between losing a district or keeping a district, it's important to compare multiple scenarios.

Currently, I am generating 10,000,000 posterior samples NMC=10000000 while only keeping 1 out of 10 THIN=10 of the samples for the output dataset OUTPOST=MCMC_EST. The first 1,000 samples are being thrown out NBI=1000. The PARMS are initial values.

The MODEL statement specifies that each population is distributed normally with mean=END and the variance=VAR_END.

This is the code:

PROC MCMC DATA=EST NBI=1000 THIN=10 NMC=10000000 SEED=17760704
OUTPOST=MCMC_EST MONITOR=(_parms_ mu) DIAGNOSTICS=NONE
PROPDIST=T;
ARRAY STATE[50];
PARMS STATE: 0;
PARMS S2 1;
PRIOR STATE: ~ NORMAL(0, VAR=S2);
PRIOR S2: ~ IGAMMA(0.00001, SCALE=100000);
mu = STATE[INDEX];
MODEL END ~ NORMAL(mu, var=END_VAR);
RUN;


I have a couple of questions regarding this approach:

1. Is this approach reasonable for my goals?
2. Is/Are there anything that I might be missing or not considering?

The data:

DATA EST;
INPUT INDEX END END_VAR @@;
CARDS;
1 15.4241325 0.0002367 2 13.5640134 0.0006933 3 15.8262773 0.0013665 4 14.9405672 0.0002576 5 17.5277513 0.0001028
6 15.5865207 0.0002908 7 15.1101493 0.0001562 8 13.8189535 0.0001989 9 16.8993923 0.0006576 10 16.2084906 0.0005570
11 14.2270479 0.0002234 12 14.3967189 0.0009461 13 16.3803101 0.0000703 14 15.7319626 0.0000578 15 14.9750436 0.0000653
16 14.9089526 0.0001230 17 15.3310353 0.0001038 18 15.3706757 0.0076526 19 14.1129090 0.0002326 20 15.6475029 0.0001247
21 15.7525401 0.0002435 22 16.1093723 0.0001705 23 15.5523285 0.0000302 24 14.9272146 0.0001557 25 15.6470303 0.0000951
26 13.8905054 0.0001742 27 14.4875277 0.0000478 28 14.9926327 0.0038334 29 14.1235860 0.0003578 30 16.0273223 0.0000774
31 14.5934536 0.0007836 32 16.8139648 0.0001331 33 16.1917873 0.0005469 34 13.5895981 0.0024562 35 16.2745662 0.0000071
36 15.2189476 0.0001872 37 15.2603779 0.0001952 38 16.3782608 0.0000426 39 13.8720971 0.0003982 40 15.4665904 0.0003467
41 13.7034732 0.0002322 42 15.7491440 0.0001953 43 17.2153554 0.0001408 44 15.0053915 0.0005280 45 13.3556029 0.0000774
46 15.9944187 0.0001168 47 15.8471508 0.0001357 48 14.4339610 0.0001527 49 15.5910714 0.0000591 50 13.3356746 0.0013045
;
RUN;

• The answer to question 2. is yes, you are missing the description of your data and of your problem. Please, give details. – utobi Dec 1 '16 at 15:02
• You cannot dump your code on readers and expect them to translate it into a meaningful model. – Xi'an Dec 1 '16 at 15:39
• As a comment: thinning is rather not useful: stats.stackexchange.com/questions/171602/… – Tim Dec 2 '16 at 12:34