I have been given an AR model with seasonal variation.  \begin{equation} (1-\theta_1B)(1-\theta_2B^8)(log(Y_t)-\mu)=\epsilon_t \end{equation} Setting $X_t=log(Y_t)-\mu$ one gets the following  \begin{equation}X_t=\epsilon_t+\theta_1X_{t-1}+\theta_2X_{t-8}-\theta_1\theta_2X_{t-9}. \end{equation}  We are given $\theta_1$, $\theta_2$ and the last 16 observations of $log(Y_t)$. $\epsilon_t$ is white noise and $\epsilon_t\sim N(0,\sigma_{\epsilon}^2)$.  The one and two step predictions thus become (where the know observations are called $\textbf{X}_t$). \begin{align} E[X_{t+1}|\textbf{X}_t]&amp;=E[\epsilon_{t+1}|\textbf{X}_t]+\theta_1E[X_{t}|\textbf{X}_t]+\theta_2E[X_{t-7}|\textbf{X}_t]-\theta_1\theta_2E[X_{t-8}|\textbf{X}_t]\\&amp;=\theta_1X_{t}+\theta_2X_{t-7}-\theta_1\theta_2X_{t-8}=\hat{X_{t+1}}\\ E[X_{t+2}|\textbf{X}_t]&amp;=E[\epsilon_{t+2}|\textbf{X}_t]+\theta_1E[X_{t+1}|\textbf{X}_t]+\theta_2E[X_{t-6}|\textbf{X}_t]-\theta_1\theta_2E[X_{t-7}|\textbf{X}_t]\\&amp;=\theta_1E[X_{t+1}|\textbf{X}_t]+\theta_2X_{t-6}-\theta_1\theta_2X_{t-7}=\hat{X_{t+2}} \end{align}  And the variance of the prediction errors become  \begin{align} V[X_{t+1}-\hat{X_{t+1}}]&amp;=V[\epsilon_{t+1}+\theta_1X_{t}+\theta_2X_{t-7}-(\theta_1\theta_2X_{t-8}-\theta_1X_{t}+\theta_2X_{t-7}-\theta_1\theta_2X_{t-8})]\\ &amp;=\sigma_\epsilon^2\\ V[X_{t+2}-\hat{X_{t+2}}]&amp;=V[\epsilon_{t+2}+\theta_1X_{t+1}+\theta_2X_{t-6}-\theta_1\theta_2X_{t-7}-(\theta_1E[X_{t+1}|\textbf{X}_t]+\theta_2X_{t-6}-\theta_1\theta_2X_{t-7})]\\ &amp;=V[\epsilon_{t+1}+\theta_1(X_{t+1}-\hat{X_{t+1}})]=(1+\theta_1^2)\sigma_\epsilon^2 \end{align}  My issue is now that I need the 95% prediction interval. What I thought (and my question is whether this is correct) I could was the following where $\hat{log(Y_{t+1})}=\hat{X_{t+1}}+\mu$, $\hat{log(Y_{t+2})}=\hat{X_{t+2}}+\mu$, $V[log(Y_{t+1})-\hat{log(Y_{t+1})}]=\sigma_\epsilon^2$ and $V[log(Y_{t+2})-\hat{log(Y_{t+2})}]=(1+\theta_1^2)\sigma_\epsilon^2$. As $log(Y_t)$ can be written as an MA process then $log(Y_t)$ follows a normal distribution and thus $Y_t$ follows a log-normal distribution. Thus using expressions from  gives that the 95% ($\alpha=0.05$) prediction interval for $Y_{t+k}$ becomes \begin{align} &amp;[lnN_{\frac{\alpha}{2}}(m_k,v_k),\ lnN_{\frac{1-\alpha}{2}}(\text{exp}(m_k,v_k)]\\ m_k&amp;=\text{exp}(\hat{log(Y_{t+k})}+\frac{V[log(Y_{t+k})-\hat{log(Y_{t+k})}]^2}{2}\\ v_k&amp;=(\text{exp}(V[log(Y_{t+k})-\hat{log(Y_{t+k})}])-1)\text{exp}(2\hat{log(Y_{t+k})}+V[log(Y_{t+k})-\hat{log(Y_{t+k})}]^2) \end{align} Where $lnN_{\frac{\alpha}{2}}(m_k,v_k)$ is the $\frac{\alpha}{2}$ quantile of a log-normal distribution with mean and variance respectively $m_k$ and $v_k$.
123722
I am attempting to use simple linear regression to construct a 95% prediction interval for a continuous response variable (Y) using a continuous input variable (X). When examining my data, I realized that if I log-transform both X and Y, the assumptions of equal variance and normality are essentially met, but if I do not log transform both variables, they are not. I'm able to run all of the R code needed to get the prediction interval for the transformed data, and I've read online and elsewhere how to interpret this interval. However, is there a way to get a 95% prediction interval on the untransformed scale using the 95% prediction interval on the log-transformed scale? In other words, I'm interested in talking about the conditional mean response instead of the median response (which, according to what I've read, is how I must interpret my prediction interval). I could exponentiate the interval endpoints, but I'm fairly sure that's incorrect because E(e^{Y}) is usually not equal to e^{E(Y)}.
74773
I have implemented Expectation maximization to fit some of the parameters of a linear Gaussian state space model using Kalman filtering / smoothing.  The model is:  $x(t) = Ax(t - 1) + w(t); w(t) \sim \mathcal{N}(0, Q)$  $y(t) = Hx(t) + v(t); v(t) \sim \mathcal{N}(0, R)$  $m(0) = m_0, P(0|0) = P_0$  // Initial mean and variance  The parameters $\theta$ consist of some of the entries of $A, Q, R$ as well as the initial conditions $m_0, P_0$.  On each iteration of the EM algorithm I (1) set the new initial conditions to $m_0 = \hat{m}(0|T), P_0 = \hat{P}(0 | T)$, where $\hat{m}, \hat{P}$ are the posterior estimates from the smoothing step.  and (2) use a numerical optimizer to maximize $\mathcal{Q}(\theta)$ (the expected log likelihood).  From what I can see on small 2D examples, this step is correctly computing the maximum of $\mathcal{Q}$.  The new initial conditions and new parameters are then fed back into the smoother to do it all again.  To check the correctness of my implementation, I have simulated data from a model with known parameters $\theta_{\text{true}}$ for a variety of sample sizes $N_1, N_2, ...$ (e.g. 100, 1000, ...) and for each of these datasets I run my EM algorithm to convergence to get parameter estimates $\hat{\theta}_1, \hat{\theta}_2, ...$  I am consistently observing that $\hat{\theta}_i \rightarrow \theta_{\text{true}}$, and moreover, this convergence is monotone.  So, it seems quite clear that my implementation is correctly inferring the parameters, and it's estimates improve with more data.  I am also looking at the log-likelihood directly, which can be calculated from the filtering recursion, and this function is increasing monotonically.  So, running my EM implementation is maximizing the likelihood itself as well.  However, during the EM iterations $t = 1, 2, ...$ the function $\mathcal{Q}(\theta^t)$ is not necessarily monotone in $t$ as it should be.  It seems to either increase monotonically (as it should), decrease monotonically, increase and then fall back down before converging, or decrease and then climb back up before converging.  In each case however, parameter estimates are still accurate.  Does anyone have suggestions for how this could be the case?  I have checked and rechecked my code and my math an innumerable number of times, and am not finding any bugs.
50986
This question regards the problem of . Vectors and matrices will be denoted in bold.  Premises. Let $N,K$ be given integers, with $K \gg N &gt; 1$. The transpose of matrix $\mathbf{A}$ will be denoted with $\mathbf{A}^T$. Suppose the following statistical model holds $$ (*) \quad \mathbf{y} = \mathbf{Hx + n}, \quad \mathbf{n} \sim \mathcal{N}_{K}(\mathbf{0}, \mathbf{C})  $$ where $\mathbf{y} \in \mathbb{R}^{K \times 1}$ are the observables, $\mathbf{H} \in \mathbb{R}^{K \times N}$ is a known full-rank matrix, $\mathbf{x} \in \mathbb{R}^{N \times 1}$ is a deterministic vector of unknown parameters (which we want to estimate) and finally $\mathbf{n} \in \mathbb{R}^{K \times 1}$ is a disturbance vector (noise) with a known (positive definite) covariance matrix $\mathbf{C} \in \mathbb{R}^{K \times K}$. The  (ML) estimate of $\mathbf{x}$, denoted with $\hat{\mathbf{x}}_{ML}$, is given by $$ (1) \quad \hat{\mathbf{x}}_{ML} = (\mathbf{H}^T \mathbf{C^{-1}} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} \mathbf{y} $$ and this is also the standard formula of Generalized Linear Least Squares (GLLS). Consider the standard formula of Ordinary Least Squares (OLS) for a linear model, i.e. $$ (2) \quad \hat{\mathbf{x}}_{OLS} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T \mathbf{y} $$ As a final note on notation, $\mathbf{I}_K$ is the $K \times K$ identity matrix and $\mathbf{O}$ is a matrix of all zeros (with appropriate dimensions). Now, for the problem at hand, assume that $\mathbf{C}^{-1} = \mathbf{I}_K + \mathbf{X}$, where $\mathbf{X} \in \mathbb{R}^{K \times K}$ is a symmetric, invertible matrix (and, for the formalism to make sense, $\mathbf{X}$ is such that $\mathbf{I}_K + \mathbf{X}$ is invertible and the inverse is positive definite). In this case, it can be proven (using matrix inversion lemma) that $$ (3) \quad (\mathbf{H}^T \mathbf{C^{-1}} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} = (\mathbf{H}^T \mathbf{H})^{-1} \mathbf{H}^T + \mathbf{Q} $$ where the expression for the matrix $\mathbf{Q} \in \mathbb{R}^{N \times K}$ can be found (I will omit it here). Now, finally,  Proposition 1. If $\mathbf{H}^T\mathbf{X} = \mathbf{O}_{N,K}$, then equation $(1)$ degenerates in equation $(2)$, i.e., there exists no difference between GLLS and OLS.  The proof is straigthforward and is valid even if $\mathbf{X}$ is singular. Now, my question is  Question: Can an equation similar to eq. $(3)$ (which "separates" an OLS-term from a second term) be written when $\mathbf{X}$ is a singular matrix?  I am not interested in a closed-form of $\mathbf{Q}$ when $\mathbf{X}$ is singular. But I do am interested in understanding the concept beyond that expression: what is the actual role of $\mathbf{Q}$? In which space does it operate?  I found this problem during a numerical implementation where both OLS and GLLS performed roughly the same (the actual model is $(*)$), and I cannot understand why OLS is not strictly sub-optimal. I found this slightly counter-intuitive, since you know a lot more in GLLS (you know $\mathbf{C}$ and make full use of it, why OLS does not), but this is somehow "useless" if some conditions are met. What are these conditions?  As a final note, I am rather new to the world of Least Squares, since I generally work within a ML-framework (or MMSE in other cases) and never studied the deep aspects of GLLS vs OLS, since, in my case, they are just intermediate steps during the derivation of MLE for a given problem.  If the question is, in your opinion, a bit too broad, or if there is something I am missing, could you please point me in the right direction by giving me references? Preferably well-known books written in standard notation.
115383
I am struggling to understand the derivation of the expected prediction error per below (ESL), especially on the derivation of 2.11 and 2.12 (conditioning, the step towards point-wise minimum). Any pointers or links much appreciated.   Below I am reporting the excerpt from ESL pg. 18. The first two equations are, in order, equation 2.11 and 2.12.    Let $X \in \mathbb{R}^p$ denote a real valued random input vector, and $Y \in \mathbb{R}$ a real valued random output variable, with joint distribution $\text{Pr}(X,Y)$. We seek a function $f(X)$ for predicting $Y$ given values of the input $X$. This theory requires a loss function $L(Y,f(X))$ for penalizing errors in prediction, and by far the most common and convenient is squared error loss: $L(Y,f(X))=(Y-f(X))^2$. This leads us to a criterion for choosing $f$,  $$ \begin{split} \text{EPE}(f) &amp;= \text{E}(Y - f(X))^2\\ &amp; = \int [y - f(x)]^2 \text{Pr}(dx, dy) \end{split} $$  the expected (squared) prediction error. By conditioning on $X$, we can write EPE as  $$ \text{EPE}(f) = \text{E}_X \text{E}_{Y|X}([Y-f(X)]^2|X) $$  and we see that it suffices to minimize EPE point-wise:  $$ f(x) = \text{argmin}_c \text{E}_{Y|X}([Y-c]^2|X) $$  The solution is  $$ f(x) = \text{E}(Y|X=x) $$  the conditional expectation, also known as the regression function.
100157
I would like to know more details about R issues mentioned in .   For e.g. the first problem still exists in R version 3.0.1  # generate an AR(1) with mean 50 set.seed(66)      # so you can reproduce these results x = arima.sim(list(order=c(1,0,0), ar=.9), n=100) + 50    mean(x)     [1] 50.60668   # the sample mean is close arima(x, order = c(1, 0, 0))     Coefficients:            ar1  intercept  &lt;--  here's the problem         0.8971    50.6304  &lt;--  or here, one of these has to change   s.e.  0.0409     0.8365   *direct copy from the webpage  Theoretical model of above simulation is: $X_t = 5 + 0.9 X_{t-1} + Z_t$, where $\{ Z_t\} $ is the white noise  Is R following particular model convention or is R wrong?  If it is wrong, it is still in R core package and I wonder why it isn't changed.  Also, it would be very helpful if I get some hint on other five issues.
41894
According to Theorem 8.8 in Time Series A.W. van der Vaart an ARMA process $$\phi (L)X_t=\theta(L)\epsilon_t$$ has a  unique stationary solution $X_t=\psi(L)\epsilon_t$ with $\psi=\theta/\phi$ if $\phi$ has no roots on the complex unit circle. This would imply that the explosive process, with $\rho&gt;1$, is a stationary process $$X_t=\rho X_{t-1}+\epsilon_t$$ with stationary solution $X_t=\sum_{i=1}^\infty \rho^{-i}\epsilon_{t+i}$.  Now indeed $\sum_{i=1}^{\infty} \rho^{-i} &lt; \infty$ so that weak stationarity can be proved by using this representation.   However, here on stackexchange I see a lot of question/answers that suggest that the process above is not stationary (see for example , ). In particular, the accepted answer of the latter question claims that the process is non-stationary by simulating a series and showing it displays explosive trending behaviour.  I think the only way to reconcile the theorem I mention above and the plots in the accepted answer of () is the following: the explosive process is indeed stationary but non-ergodic, that is, we cannot find the statistical properties of $X_t$ such as $\mathbb{E}(X_t)=\mu$ by observing a single infinitely long sample path of the explosive process, mathematically: $$\lim_{t \to \infty}\frac{1}{t}\sum_{t=1}X_t \neq\mathbb{E}X_t$$  Is this reading correct?
53362
I am a self learner, and I am studying time series analysis.  I came through the fact that ARIMA  can be used  to model a time series which is not stationary (Integrated ARMA model). The non stationarity taken in to account in an ARIMA model is the nonstationarity that can be removed by differencing.   I am wondering whether a  seasonality is the type of nonstationarity that can be removed by differencing?    If yes, so for a non seasonal data, we may use ARIMA modeling, so why we have SARIMA modeling ?   If no, how seasonality can be  removed in order to make data stationary?
10661
I am trying to forecast time series data for next $n$ time periods and I am using various forecasting techniques like simple moving average, exponential smoothing (Single,double,triple) and auto.arima().  I am getting equal forecasting values for simple moving average and single exponential smoothing and at times, even auto.arima() function in R is giving me the same forecast values.  Here is the code that I am trying to use  mod_mva &lt;- SMA(data.train,n=2) forecast_mva &lt;-forecast(mod_mva,ts_freq)   Here is the output:
64733
I have created a bayesian model that estimates 6 parameters using rjags from R. Now i want to do some predictions based on new data in R. Can anyone help me with an example.  model {    B.m[1] &lt;- Kb # mean B[1]    B[1] ~ dnorm( B.m[1], tau_B)    for(t in 2 : 10) { # Process model       B.m[t]&lt;- max( B[t-1] + rB[t-1](1 -B[t-1]/K) - C[t-1] , 0.0)       B[t] ~ dnorm(B.m[t], tau_B) } # Observation mode       for (t in 1: 10) {          LogC_mean[t] &lt;- log( max(qsegB[t]Eff[t], 0.000001))          C[t] ~ dlnorm(LogC_mean[t], tau_c) } # Prior          b~dbeta(1, 1)          cv &lt;- 1          K.mean &lt;- 500000          K.v &lt;- pow(cvK.mean, 2)          K.tau &lt;- 1/K.v          K ~ dnorm(K.mean, K.tau)          r ~ dlnorm(0.0, 1.4)          qseg ~ dlnorm(-15.4, 1.44)          tau_c ~ dgamma(0.001,0.001)          tau_B ~ dgamma(0.01,0.01)        }     }
58121
I'm using R(3.1.1), and ARIMA models for forecasting. I would like to know what should be the &quot;frequency&quot; parameter, which is assigned in the ts() function, if im using time series data which is:  separated by minutes and is spread over 180 days (1440 minutes/day) separated by seconds and is spread over 180 days (86,400 seconds/day).  If I recall right the definition, a &quot;frequency&quot; in ts in R, is the number of observations per &quot;season&quot;. Question part 1: What is the &quot;season&quot; in my case? If the season is &quot;day&quot;, then is the &quot;frequency&quot; for minutes = 1440 and 86,400 for seconds? Question part 2: Could the &quot;frequency&quot;  also depend on what I am trying to achieve/forecast? for example, in my case, I'd like to have a very short-term forecast. One-step ahead of 10minutes each time. Would it then be possible to consider the season as an hour instead of a day? In that case frequency= 60 for minutes, frequency = 3600 for seconds? I've tried for example to use frequency = 60 for the minute data and got better results compared to frequency = 1440 (used fourier see link below by Hyndman)  (The comparison was made by using MAPE for the measure of forecast accuracy) In case the results are complete arbitrary, and the frequency cannot be changed. What would be actually the interpretation of using freq = 60 on my data? I also think it's worth mentioning that my data contains seasonality at every hour and every two hours (by observing the raw data and the Autocorrelation function)
2768
I have scoured around, reading posts on Cross Validated () and also looking at references including Dobson and McCullagh and Nelder, e.g.  so I am aware that this topic is well trodden. Nevertheless, I am trying to articulate and formalize my understanding of GLM and while several posts have helped me with that, I am conscious of gaps and the possibility of an unsound foundation of my understanding.   In simple linear regression we have some set of observations $(x_i, y_i)$ pairs and treat $y_i$ as a realization of a random variable, $Y_i$, distributed as $Y_i \sim N(\mu_i, \sigma^2)$. The means ($\mu_i$) depend on the predictor but the variance is constant. We model $\mu_i = \beta_0 + (\beta_1 \times x_i)$ (or matrix equivalent ) which I believe is the same as saying $Y_i = \mu_i + \epsilon_i$ with $\epsilon_i \sim N(0,\sigma^2)$. I am not sure but I think a correct way to state $\mu$ is $\mu_i = E[Y_i|X=x]$. Can someone confirm that?  Anyway, we might transform the response to achieve linearity (maybe taking logs) and in that case we are modelling $log(Y_i) = \alpha_0 + (\alpha_1 \times x_i) + \epsilon_i$ where $\epsilon_i$ now is assumed to have a log normal distribution.   We generalize by decomposing the model into:   A structural component ($\beta_0 + (\beta_1 \times x_i)$) A link g(.) A response distribution (or random component) (member of the exponential family - Guassian, bin,gamma etc)   Let's say that the observations are assumed to come from a distribution in the exponential family and to keep things simple assume that it is the Gaussian distribution. Again, the expected value comes out as $E[Y_i|X=x] = \mu_i$, but does so from the first derivative of the b(theta) term in the exponential family form () (pg 29). In the case where we assume a binomial distribution, the expected value comes out at np but it still comes from the first derivative of the b(theta) term when the distribution is expressed in the exponential family form. Note, I don't fully appreciate or understand the derivation of the mean or variance at this stage, perhaps someone could provide a layman's explanation?  Instead of modelling the mean as was done for simple linear regression, we now model a transformation of the mean so instead of saying $\mu_i = \beta_0 + (\beta_1 \times x_i)$ we are saying $g(\mu_i) = \eta_i = \beta_0 + (\beta_1 \times x_i)$ where g is some link function (invertible and differentiable). I think this is a key (yet still somewhat confusing to me) distinction between SLR and GLM. In SLR we transform the response ($y_i$) and model that, in GLM we transform the expected value ($\mu_i$ in the case of the Gaussian example) and model that. Another way of saying it is that in the SLR case we model $E[g(Y_i|X=x)] = \beta_0 + (\beta_1 \times x_i)$ but in the GLM world we are modelling $\eta_i = g[E(Y_i|X=x)] = \beta_0 + (\beta_1 \times x_i)$.   My question is around verifying my understanding and statement of the essence of the foundations of GLM and the differences between GLM and the traditional linear model. Thanks.
13876
I would appreciate if someone could help me write the mathematical equation for the seasonal ARIMA (2,1,0) x (0,2,2) period 12. I'm a little confused with how to go about this. I would prefer an equation involving $Y_t , e_t, θ$ and $Θ$.
56710
Given a random walk $x_t$,   $$x_t=x_{t-1}+\varepsilon_t,$$  consider estimating the slope coefficient $\beta$ in   $$x_t=\beta x_{t-1}+\varepsilon_t$$  by OLS.  question and the following answer noted that $\hat \beta^{OLS}$ has a skewed distribution. My question is,  Are any of OLS assumptions violated in the model $x_t=\beta x_{t-1}+\varepsilon_t$ given that $x_t$ is a random walk? If so, what are the violations?
25216
suppose for a time series, the auto.arima output is:  ARIMA(1,0,1) with non-zero mean   Coefficients:          ar1     ma1  intercept       0.7456  0.2775   767.7463 s.e.  0.0804  0.1197     0.1072   It is specified in help file of arima function that definition used for this function is  X[t] = a[1]X[t-1] + … + a[p]X[t-p] + e[t] + b[1]e[t-1] + … + b[q]e[t-q]   and if we include mean(m),  this formula applies to X - m rather than X.  My questions are:  (1)  For the model considered above, what will be the arima equation?    (a) X[t] =  767.7463 + 0.7456*X[t-1] + 0.2775*e[t-1]    or    (b) X[t] =  195.3258 + 0.7456*X[t-1] + 0.2775*e[t-1],   where constant term is achieved by 195.3258 = 767.7463*(1-0.7456).  (2) How to interpret this constant term if we have more than one AR coefficient?  e.g. ARIMA(2,0,2) with non-zero mean   Coefficients:           ar1      ar2     ma1     ma2  intercept       -0.0124  -0.6945  0.3629  0.9414    -0.1288 s.e.   0.1172   0.1132  0.0598  0.1270     0.0636   (3) What should be the ARIMA equation in presence of drift term?  e.g. ARIMA(1,1,0) with drift           Coefficients:          ar1   drift       0.3307  0.3981 s.e.  0.0946  0.1463
495
This is a basic question on Box-Jenkins MA models. As I understand, an MA model is basically a linear regression of time-series values $Y$ against previous error terms $e_t,..., e_{t-n}$. That is, the observation $Y$ is first regressed against its previous values $Y_{t-1}, ..., Y_{t-n}$ and then one or more $Y - \hat{Y}$ values are used as the error terms for the MA model.  But how are the error terms calculated in an ARIMA(0, 0, 2) model? If the MA model is used without an autoregressive part and thus no estimated value, how can I possibly have an error term?
112658
Is it enough to show that MSE = 0 as $n\rightarrow\infty$? I also read in my notes something about plim. How do I find plim and use it to show that the estimator is consistent?
226
On page 19 of the textbook  (by James, Witten, Hastie and Tibshirani--it is freely downloadable on the web, and very good), the following is stated:  Consider a given estimate $$\hat{Y} = \hat{f}(x)$$ Assume for a moment that both $$\hat{f}, X$$ are fixed. Then, it is easy to show that: $$\mathrm{E}(Y - \hat{Y})^2 = \mathrm{E}[f(X) + \epsilon - \hat{f}(X)]^2$$ $$ = [f(X) - \hat{f}(X)]^2 + \mathrm{Var}(\epsilon)$$  It is further explained that the first term represents the reducible error, and the second term represents the irreducible error. I am not fully understanding how the authors arrive at this answer. I worked through the calculations as follows: $$\mathrm{E}(Y - \hat{Y})^2 = \mathrm{E}[f(X) + \epsilon - \hat{f}(X)]^2$$ This simplifies to $[f(X) - \hat{f}(X) + \mathrm{E}[\epsilon]]^2 = [f(X) - \hat{f}(X)]^2$ assuming that $\mathrm{E}[\epsilon] = 0$. Where is the $\mathrm{Var}(x)$ indicated in the text coming from? Any suggestions would be greatly appreciated.