Omitted Variable Bias, verification in Gretl I am trying to verify the expression for Omitted Variable Bias (OVB) as given e.g. in Wooldridge: $\tilde{\beta_1} = \hat{\beta_1} + \hat{\beta_2} \cdot \tilde{\delta_1}$, where $\tilde{\delta_1}$ is the estimated slope of the regression of $x_2$ on $x_1$.
Choosing the housing price data (hprice1.gdt) from Wooldridge available for gretl, I obtain the following estimates for the relevant regression coefficients:
Model1 (price ~ sqft)
  ------------------------
  const       11.2041         
  sqrft       0.140211                

Model2 (price ~ sqft + bdrms)
  ------------------------
  const       -19.315       
  sqrft       0.128436                
  bdrms       15.1982                 

Model3 (bdrms ~ sqft)
  ------------------------
  const       2.00808       
  sqrft       0.000774748             

So $\tilde{\beta_1}=0.140211$, $\tilde{\delta_1} = 0.000774748$, $\hat{\beta_2}=15.1982$ and $\hat{\beta_1}=0.128436$
But these values do not quite match the expression from above:
$0.124836+15.1982*0.000774748 = 0.1366108 \neq 0.140211$
Where is my misunderstanding ? Is the formula above valid only in the limit of infinite sample size?
 A: As pointed out in the comments to an older incorrect version of this answer, the relation is exact.
We assume the correct specification is (matrix notation, lowercase is column vector)
$$ y = X_1\beta_1 + \beta_2x_2 + \varepsilon \tag{1}$$
where we have subsumed the series of $1$'s in $X_1$ (which is now a matrix with two columns, while $\beta_1$ is a column vector with two columns including also the constant), and with the error term being mean-independent from the regressors, $E(\varepsilon\mid X_1, x_2) = 0$.
Assume that we run this regression. We will obtain coefficient estimates and a series of residuals. Then, copying from the comment, it holds that 
$$y = X_1 \hat{\beta}_1 + x_2 \hat{\beta}_2  + \hat{\varepsilon}  \tag{2}$$
Assume now that we ignore the $x_2$ regressor and we estimate
$$y = X_1\beta_1 + u \tag{3}$$
then $\tilde \beta_1 = \left(X_1'X_1\right)^{-1}X_1' y $ and using $(2)$
$$\tilde \beta_1 = \left(X_1'X_1\right)^{-1}X_1' X_1 \hat{\beta}_1 + \left(X_1'X_1\right)^{-1}X_1'x_2 \hat{\beta}_2 + \left(X_1'X_1\right)^{-1}X_1'\hat{e} \tag{4}$$
But by construction, $\left(X_1'X_1\right)^{-1}X_1'\hat{e} = 0$, so we are left with 
$$\tilde \beta_1 = \hat{\beta}_1 + \left(X_1'X_1\right)^{-1}X_1'x_2 \hat{\beta}_2 \tag{5}$$
Then note that $\left(X_1'X_1\right)^{-1}X_1'x_2 = \tilde \delta_1$ from the specification 
$$x_2 = X_1\delta_1 + \epsilon \tag{6}$$
so we end up with 
(for which we assume that there are no other variables that covary with both $X_1$ and $x_2$ and so that $\tilde \delta_1$ is a consistent estimator), and so
$$\tilde \beta_1 = \hat \beta_1 + \tilde \delta_1 \hat \beta_2 \tag{7}$$
This is exact. Since in all the above calculations $\beta_1$ and $\delta_1$ include two coefficients, the exact relationship is 
$$(\tilde \beta_{10}, \tilde \beta_{11}) = (\hat \beta_{10}, \hat \beta_{11}) + (\tilde \delta_{10}, \tilde \delta_{11})\hat \beta_2$$ 
where the second zero subscript indicates the estimate for the constant term.
So we have to check 
$$\tilde \beta_{10} = 11.2041 = \hat \beta_{10} + \tilde \delta_{10}\hat \beta_2 = -19.315 + 2.00808\cdot 15.1982 = 11.2042$$
$$\tilde \beta_{11} = 0.140211 = \hat \beta_{11} + \tilde \delta_{11}\hat \beta_2 = 0.128436 + 0.000774748\cdot 15.1982 = 0.1402107$$
and we're good. The OP made a simple mistake in punching the numbers in the calculations (he used $0.124836$ instead of the correct one $0.128436$).
A: 
Where is my misunderstanding ? Is the formula above valid only in the
  limit of infinite sample size?

The formula is valid and exact for all sample sizes. Your misunderstanding is a simple typo, there's no misunderstanding at all.
When  you wrote:

But these values do not quite match the expression from above:
  $0.124836+15.1982*0.000774748 = 0.1366108 \neq 0.140211$

Somehow you put $0.124836$ instead of $0.128436$ switching the $4$ for the $8$. Fixing the typo  gives you the expected result:
$$0.128436+15.1982*0.000774748=0.140211$$
The proof is rather simple. Let $Y$ denote price,  $X$ denote sqrft and $Z$ denote bdrms. Then:
$$
\tilde{\beta}_1 = \frac{cov(X, Y)}{var(X)}= \frac{cov(X, \hat{\beta}_1X + \hat{\beta}_2Z + \hat{\epsilon})}{var(X)} = \hat{\beta}_1 + \hat{\beta}_2\frac{cov(X, Z)}{var(X)} = \hat{\beta}_1+ \hat{\beta}_2\tilde{\delta}_1
$$
Where we know $cov(X, \hat{\epsilon}) = 0$ by construction in OLS and $\frac{cov(X, Z)}{var(X)}$ is the coefficient of regressing $Z \sim X$ which we are denoting for $\tilde{\delta}_1$.
This relationship is exact and just a simple property of the algebra of OLS.
If you want to manually check this in R, there's a package called wooldridge with all the datasets from the textbook:
library(wooldridge)
data("hprice1")
coef(lm(price ~ sqrft, hprice1))[2]
#  sqrft 
# 0.140211 
coef(lm(price ~ sqrft + bdrms, hprice1))[2] + 
  coef(lm(price ~ sqrft + bdrms, hprice1))[3]*coef(lm(bdrms ~ sqrft , hprice1))[2]
#  sqrft 
# 0.140211 

A: I'm looking at the Wooldridge textbook 2006 international version, where on p.96 in equ. 3.45 it says: $E(\tilde{\beta_1}) = E(\hat{\beta_1})+E(\hat{\beta_2}) \tilde{\delta_1}$.
This is compatible with the previous answer by Alecos Papadopoulos, but it is not (in general) the equation that you stated in your question, without expectation operators. That is, this relation only holds in expectation, not in any given sample. You would have to conduct a Monte Carlo simulation exercise to calculate the expectations involved (up to arbitrary precision by increasing the number of draws) and then the stated equality should show up pretty much exactly.
