Short answers:
1. It simplifies it. (Frankly, I did not get the question).
2. No, you can never ignore it, as lack of i.i.d. has immediate consequences on the variances of whatever you are estimating.
Medium answer: Pretty much the central issue with the bootstrap is, 'Does the proposed procedure reproduce the features of the data?'. Violation of the i.i.d. assumption is a big deal: your data are dependent, you (most likely) have less information in your data than you would have in an i.i.d. sample of the same size, and if you run a naive bootstrap (resample the individual observations), the standard errors you get from it will be too small. The proposed procedure circumvents the problem of lack of independence by capturing (or at least attempting to capture) the dependence in the model structure and parameters. If successful, each bootstrap sample would reproduce the features of the data, as needed.
Long answer: There are multiple layers of assumptions concerning the bootstrap, and even in the simplest possible case (i.i.d. data, estimation of the mean), you have to make at least three assumptions: (1) the statistic of interest is a smooth function of the data (true in the case of the mean, not so true even in the case of percentiles, totally off with say nearest neighbor matching estimators); (2) the distribution from which you bootstrap is "close" to the population distribution (works OK in the case of i.i.d. data; may not work OK in the case of dependent data, where you essentially have only one trajectory = one observation in the case of time series, and you have to invoke additional assumptions like stationarity and mixing to strecth this single observation into a quasi-population); (3) your Monte Carlo bootstrap sampling is a good enough approximation to the complete bootstrap with all possible subsamples (the inaccuracy from using Monte Carlo vs. the complete bootstrap is much less than the uncertainty you are trying to capture). In the case of the parametric bootstrap, you also make an assumption that (4) your model perfectly explains all the features of the data.
As a warning of what could go wrong with (4), think about regression with heteroskedastic errors: $y=x\beta + \epsilon$, Var$[\epsilon] = \exp[ x\gamma]$, say. If you fit an OLS model and resample the residuals as if they were i.i.d., you will get a wrong answer (some sort of $\bar\sigma^2 (X'X)^{-1}$ where $\bar\sigma^2$ is the average $1/n \sum_i \exp[x_i \gamma]$, instead of the appropriate $(X'X)^{-1} \sum \exp[x_i \gamma] x_i x_i' (X'X)^{-1}$). So if you wanted to have a fully parametric bootstrap solution, you would've have to fit the model for heteroskedasticity along with the model for the mean. And if you suspect serial or other sort of correlation, you would've have to fit the model for that, too. (See, the non-parametric distribution-free flavor of the bootstrap is pretty much gone for now, as you have replaced the voice of the data with the synthesized voice of your model.)
The method you described works around the i.i.d. assumption by creating a whole new sample. The greatest problem with the dependent data bootstrap is to create the sample that would have the dependence patterns that would be sufficiently close to those in the original data. With time series, you could use block bootstraps; with clustered data, you bootstrap the whole clusters; with heteroskedastic regression, you have to with wild bootstraps (which is a better idea than the bootstrap of residuals, even if you have fitted a heteroskedasticty model to it). In the block bootstrap, you have to make an educated guess (or, in other words, have good reasons to believe) that distant parts of time series are approximately independent, so that all of the correlation structure is captured by the adjacent 5 or 10 observtations that form the block. So instead of resampling observations one by one, which totally ignores the correlation structure of the time-series, you resampling them in blocks, hoping that this would respect the correlation structure. The parametric bootstrap you referred to says: "Rather than fiddling with the data and assembling the new dolls from the pieces of the old ones, why don't I just stamp the whole molded Barbie for you instead? I've figured out what kind of Barbies you like, and I promise I will make you one you'd like, too."
In case of the parametric bootstrap you described, you have to be pretty damn sure that your HMM model fit is pretty much perfect, otherwise your parametric bootstrap may lead to incorrect results (Barbies that cannot move their arms). Think about the above heteroskedastic regession example; or think about fitting an AR(1) model to AR(5) data: whatever you do with the parametrically simulated data, they won't have the structure the original data used to have.
Edit: as Sadeghd clarified his question, I can respond to that, as well. There is a humongous variety of the bootstrap procedures, each addressing the particular quirk in either the statistic, the sample size, the dependence, or whatever an issue with the bootstrap could be. There is no single way to address dependence, for instance. (I've worked with survey bootstraps, there are about 8 different procedures, although some are mostly of methodological rather than practical interest; and some are clearly inferior in that they are only applicable in special, not easily generalizable, cases.) For a general discussion of issues you could face with the bootstrap, see Canty, Davison, Hinkley and Ventura (2006). Bootstrap diagnostics and remedies. The Canadian Journal of Statistics, 34 (1), 5-27.
This doesn't exactly answer your question, but it might still solve your problem of needing to calculate risk ratios. The epiR package allows you to calculate risk ratios.
I could not get your example to work (see my comment to your question), so here is an example from the package's documentation:
library(epiR) # Used for Risk ratio
library(MASS) # Used for data
dat1 <- birthwt; head(dat1)
## Generate a table of cell frequencies. First set the levels of the outcome
## and the exposure so the frequencies in the 2 by 2 table come out in the
## conventional format:
dat1$low <- factor(dat1$low, levels = c(1,0))
dat1$smoke <- factor(dat1$smoke, levels = c(1,0))
dat1$race <- factor(dat1$race, levels = c(1,2,3))
## Generate the 2 by 2 table. Exposure (rows) = smoke. Outcome (columns) = low.
tab1 <- table(dat1$smoke, dat1$low, dnn = c("Smoke", "Low BW"))
print(tab1)
## Compute the incidence risk ratio and other measures of association:
epi.2by2(dat = tab1, method = "cohort.count",
conf.level = 0.95, units = 100, outcome = "as.columns")
Best Answer
I'm not sure what your new data.frame looks like, but you need to make sure they have the same variable names as the variables in your model. So in your model, your variable names are y and x. Then
new = data.frame(X,Y)
colnames(new) = c("x","y")
nb_pred <- predict(nb_model, new_data=new, type="response", se.fit=T)