Combining sources of uncertainty/variation in a multi-layered linear model I’m having trouble understanding how to combine parameter uncertainty and interannual variability from various levels in a bootstrapped linear model. Specifically, this model is designed to generate estimates of some phenomenon $Y$ for 0.5º grid cells in a spatial region (e.g, Central America).
It looks like this:
$$
    \boldsymbol{Y_{m,y}} = \boldsymbol{\beta_{m,y}X_{m,y}}
$$
where $\boldsymbol{\beta}$ is the vector of parameters being estimated (size 1x3), $\boldsymbol{X}$ is a matrix of observations of land cover fractions (size 3xN), $\boldsymbol{Y}$ is the vector of predicted values (size 1xN), and the subscripts indicate that the model is being parameterized for each month $m$ and year $y$ (e.g., May 2006, June 2006, and May 2007 can all have different parameter estimates). (No intercept, although I don’t think that matters for my question.)
If we think of the model as $Y_{m,y}=\sum\limits_{i=1}^3{\beta_{i,m,y}X_{i,m,y}}$, then each term $\beta_iX_i$ can be considered a subset of $Y$ — i.e., 
$Y_i$, the amount of $Y$ that happens on land cover type $i$  — in a real physical sense. Thus, there will likely be some correlation among the various parameter estimates. This is embodied/exacerbated by the fact that for each grid cell $g$, $\sum\limits_{i=1}^3{\boldsymbol{X_{g,i,m,y}}} = 1$ — that is, the land cover type fractions add up to 1.
To create uncertainty bounds for the parameters, I have done 10,000 bootstrapping runs. Let’s say there are 500 grid cells in the region. For the first bootstrapping run, I randomly choose 500 grid cells from that region with replacement. I then estimated the parameters for each month and year using that sample set. I then repeated this sampling-fitting procedure 9,999 more times.
My goal is to compare the measured/observed mean annual amount of phenomenon $Y$ with the estimated/modeled amount. I would like the uncertainty bars around the measured amount (which will take into account just interannual variability in the measured amount) to overlap with the uncertainty bars around the estimated amount (which will take into account both uncertainty in the model parameters as well as interannual variability in the annual estimated amount).
For the observations, it’s easy enough — sum up observed $Y$ for all the months in each year, then find the standard deviation across all the years. I get confused when thinking about how to do this for the estimates, though. The problem becomes even stickier if I want to compare observed vs. estimated annual GLOBAL $Y$ (i.e., across all regions).
I know that adding standard deviations in quadrature is going to be key here, but the fact that there are so many levels is confusing me. I'm working in Matlab now and am also fluent in R, but any help you can provide would be much appreciated.
 A: I have an idea for how to do this (this is Idea 1; see Idea 2 in a different answer). I won't mark this as the correct answer until I get some confirmation, so please let me know if this looks right to you (although make sure you understand what I'm trying to do!). My main point of concern is the way I include covariances in quadrature -- I'm afraid I am being redundant, and therefore the estimates of SD will be too large.
Note:


*

*$cov_b(i,j)$ means "the covariance of $i$ and $j$ across all bootstrapping runs $b$"


Find the mean and SD for a given cover ($c$) in a given month ($m$) of a given year ($y$) in a given region ($r$) across all B bootstrapping runs (e.g., forest in May 2006 in Central America):
$$
\overline{Y_{c,m,y,r}} = mean_{b}(Y_{c,m,y,r,b}) = \frac{\sum\limits_{b=1}^B{Y_{c,m,y,r,b}}}{B}
$$
$$
\sigma_{c,m,y,r} = sd_{b}(Y_{c,m,y,r,b}) = \sqrt{\frac{\sum\limits_{b=1}^B{(Y_{c,m,y,r,b}-\overline{Y_{c,m,y,r}})^2}}{B}}
$$
Find the mean and SD of total (i.e., across all 3 covers $c$) Y for a given month ($m$) of a given year ($y$) in a given region ($r$) (e.g., May 2006 in Central America):
$$
\overline{Y_{m,y,r}} = \sum\limits_{c=1}^3{\overline{Y_{c,m,y,r}}}
$$
$$
Y_{m,y,r,b} = \sum\limits_{c=1}^3{Y_{c,m,y,r,b}}
$$
$$
\sigma_{m,y,r} = \sqrt{\sum\limits_{c=1}^3{\sigma^2_{c,m,y,r}}+2\sum\limits_{c1=1}^2{\sum\limits_{c2=2\neq{c1}}^3{cov_b(Y_{c1,m,y,r,b},Y_{c2,m,y,r,b})}}}
$$
Find the mean and SD of total Y for a given month ($m$) in a given region ($r$) (e.g., for all Mays in Central America).
$$
\overline{Y_{m,r}} = mean_y(\overline{Y_{m,y,r}}) = \frac{\sum\limits_{y=1}^N{\overline{Y_{m,y,r}}}}{N}
$$
$$
Y_{m,r,b} = \frac{\sum\limits_{y=1}^N{\sum\limits_{c=1}^3{Y_{c,m,y,r,b}}}}{N}
$$
$$
\sigma_{m,r} = \frac{1}{N}\sqrt{\sum\limits_{y=1}^N{\sigma^2_{m,y,r}}+2\sum\limits_{y1=1}^{N-1}{\sum\limits_{y2=2\neq{y1}}^N{cov_b(Y_{m,y1,r,b},Y_{m,y2,r,b})}}}
$$
Find the mean and SD of annual total Y for a given region ($r$) (e.g., for Central America):
$$
\overline{Y_{r}} = \sum\limits_{m=1}^{12}{\overline{Y_{m,r}}}
$$
$$
\sigma_r = \sqrt{\sum\limits_{m=1}^{12}\sigma^2_{m,r}+2\sum\limits_{m1=1}^{11}{\sum\limits_{m2=2\neq{m1}}^{12}{cov_b(Y_{m1,r,b},Y_{m2,r,b})}}}
$$
A: I am not sure I understood the problem, so I'll star by trying to resume all important features for the discussion. Correct me if wrong, I'll edit the answer (or someone will) and we will work on it iteratively.
You have a linear model that predicts the occurrence some phenomenon Y of which you want to estimate variability. This phenomenon is calculated for a month $m$ and year $y$ and the input parameters of it are just 3: the fraction of terrain of each of the three types covered by the geographic grid. These fractions add up to one, naturally. Your objective is to estimate the inter-annual variability as measured, as estimated by the model, and to compare these with the variability that is intrinsically associated to the model.
The first thing I find hard to understand is the concept of "region". When you draw from a region, my understanding is that you are drawing groups of $\beta_i$ such that these are valid ($\Sigma_i \beta_i=1$) and that somehow are still equivalent (maybe in the error bars) of the original observed grid point. I do not understand if $Y_{m,y}$ only makes sense for the context of a wide-scale value calculated for all the grid points, resulting from the N observations on N grid points, or if it actually that runs on the gridpoints as well. In other words, I do not know what N refers to and what are the dependent variables of your model.
I think your method of calculating variability through bootstrap is justified, and you could go even for more "raw" methods, more conservative but that make less assumptions, like the jacknife, but I think before you should try to present your problem in a clearer way.  
A: I've come up with an alternative that is much simpler (this is Idea 2; see Idea 1 in a different answer). I still don't know if this is right, though! This explanation is for a given region. I will be using Matlab code to show what I'm doing.
Definitions:


*

*Z = number of years modeled

*B = number of bootstrapping runs

*estY_zbmc = $Z \times B \times 12 \times 3$ array, with model-estimated $Y$ on each land cover type $c$ in each month $m$ (e.g., July or August) in each year $z$, for each bootstrapping run $b$.

*obsY_z = $Z \times 1$ array of total $Y$ (i.e., for all land cover types combined) for each year.


Here's the code:
estY_zb = sum(sum(estY_zbmc,4),3) ; % Find total estimated Y for each year in each bootstrapping run
estY_z_mean = mean(estY_zb,2) ; % Across all bootstrapping runs, find mean total estimated Y for each year
estY_z_sd = std(estY_zb,2) ; % Across all bootstrapping runs, find st. dev. of total estimated Y for each year
estY_meanAnn = mean(estY_z_mean) ; % Find mean annual total estimated Y

That gives us the mean annual Y easily enough. What we want next is to find the error bars around that number. We'll have to do something like this:
$
\sigma = \frac{1}{Z} \sqrt{\sum\limits_{z=1}^Y{\sigma_z}+2\sum\limits_{z1=1}^{Z-1}{\sum\limits_{z2=2 \neq z1}^Z{\sigma_{z1,z2}}}}
$
The expression under the radical is just the sum of all the elements in the covariance matrix, so I can do the following (transp used because cov(X) treats rows of matrix X as observations, and columns as variables):
cov_z1z2 = cov(transp(estY_zb)) ; % transp() used because cov(X) treats rows of matrix X as observations, and columns as variables
sum_cov_z1z2 = sum(cov_z1z2(:)) ; 

And now:
estY_sd = sum_cov_z1z2 / Y ;

I can then compare estY_meanAnn and estY_sd with the mean and SD of mean annual observed $Y$:
obsY_meanAnn = mean(obsY_z) ;
obsY_sd = std(obsY_z) ;

