A simple Parametric Bootstrap on a Linear Model My model: $y_t=\beta_1+\beta_2x_{2t}+\beta_3x_{3t}+u_t$, with $u_t\sim N(0,\sigma^2)$.
My sample is of size 50.
I'm trying to do a parametric bootstrap for the t-statistic of $\beta_3$, and
I'm certain that the t-statistic is 2.349, and the respective p-value (for 47 df) is 0.0231
Since the t-distribution is symmetric around zero, i.e. $F(-x)=1-F(x)$, there are several ways to compute an approximate p-value. The one I'm using is
$$\hat p(t_{\beta_3})=2\times\frac{\sum_{i=1}^B I(t_{i}^*<-|t_{\beta_3}|)}{B}$$ where $B$ is number of bootstrap samples I generate under OLS estimates for the coeeficients of the model above, with an estimate of the variance calculated from $SSR/47$.
Below you'll find the program I'm trying to create to obtain an approximate p-value of for the t-statistic
bootstrap[ts_, B_] := Module[{counter, tsboot, ctsboot},
    counter = 1; ctsboot = 0;
    beta = lm["BestFitParameters"];
    s2 = lm["EstimatedVariance"];
    While[counter <= B,
      resboot = RandomVariate[NormalDistribution[0, s2], 50];
      yboot = Transpose[ToExpression[{x1, x2, x3}]].beta + resboot;
      lmboot = LinearModelFit[Transpose[ToExpression[{x2, x3, yboot}]],{w2,w3}, {w2, w3}];
      tsboot = lmboot["ParameterTStatistics"][[3]];
   (*Print["tsboot = ", tsboot, " and ts = ",ts];*)
      If[tsboot <- Abs[ts],(*If I remove the absolute values, I get similar results*)
      ctsboot++;
    ];
   counter++;
   ];
  ctsboot/B
  ]

To call the programme I do bootstrap[lm["ParameterTStatistics"][[3]], B], where lm["ParameterTStatistics"][[3]]=2.349 which is $t_{\beta_3}$, and $B \in \{99,999,9999\}$
The problem with my program is that the p-value computed is $0$, even for $B=9999$. Searching CV, I thought I might be creating incorrectly the bootstrap samples, but I don't see how...
I've confirmed that 'yboot', the bootstrap sample, is always different with each iteration. If for the formula of the p-value, I use $2\times\frac{\sum_{i=1}^B I(t_{i}^*>|t_{\beta_3}|)}{B}$, I get a p-value of around 0.6, for several different $B$ values
Any help would be appreciated.
 A: Your formulas are "almost" correct. First of all, if you are interested in the 2-sided hypothesis, I would recommend the following formula (rather than the two formulas you proposed):
$$ \hat{p}_{boot} = \frac{\sum_{i=1}^B I(|t_{i}^*|>|t_{\beta_3}|)}{B}$$
Now, why did I say that your formulas are "almost" correct? They would be correct, if you were correctly constructing your bootstrap estimated $t_{i}^*$, but your resampling is incorrect!
Bootstrap hypothesis testing is a bit more nuanced than bootstrapping for confidence intervals, standard errors, etc. In particular, you need to simulate from your null hypothesis [1]. In your case, this means that you should be simulating from:
$$y_t=\beta_1^*+\beta_2^* x_{2t}+ u_t$$
rather than from:
$$y_t=\beta_1^*+\beta_2^* x_{2t}+ \beta_3^* x_{3t} + u_t$$
 [1] Alternatively, you could continue simulating as you do, but then you would need to calculate p-values by inverting the confidence intervals.  
A: Based on the answer of 'air' user, here's the code that seems to be working.
bootstrap[ts_, B_] := Module[{counter, tsboot, ctsboot},
  counter = 1; ctsboot = 0;
  x1 = ToExpression[x1];
  x2 = ToExpression[x2];
  x3 = ToExpression[x3];
  y = ToExpression[y];

  (*under the Null the model has beta_ 3 = 0 *)
  X2 = Transpose[{x1, x2}];
  lmnull = LinearModelFit[Transpose[{x2, y}], {w2}, {w2}];
  beta = lmnull["BestFitParameters"][[1 ;; 2]];
  s2 = lmnull["EstimatedVariance"];

  While[counter <= B,
   resboot = RandomVariate[NormalDistribution[0, s2], 50];
   yboot = X2.beta + resboot;

   (*The bootstrap statistic is calculated for the unrestricted model \
but with yboot*)
   lmboot = 
    LinearModelFit[Transpose[{x2, x3, yboot}], {w2, w3}, {w2, w3}];
   tsboot = lmboot["ParameterTStatistics"][[3]];

   If[tsboot > Abs[ts],
    ctsboot++;
    ];
   counter++;
   ];
  2*ctsboot/B
  ]

