Which standard deviation to use for p-value of regression coefficients? Matlab generates regression coefficients for vector time series models ("vgxvarx").
For a given regression coefficient, the p-value is just the area under the t-distribution beyond |t| > (t-statistic), which is easy to do, as per this answer on Matlab Central. The expression is
pValue = 2 * tcdf(-abs(t), dfe).

I have two questions about the conversion from t to p.
The first has to do with the fact that the standard errors are described as maximum likelihood and not corrected to reflect the degrees of freedom. From Wikipedia, it seems that I would need to modify t by a factor sqrt(df/n), where n is the number of observations and df is n less the number of estimated parameters.  Can anyone chime in with how advisable this is?  My gut says that it boils down to whether I want an estimate of dispersion that is maximum-likelihood or reduced-bias.  So I fear that there is probably no right answer, and I'm hoping that someone can dispel that fear.
The second question has to do with the factor of 2 in the above expression for pValue.  I think I want the factor of 2 because I'm interested in the area under the t-distribution for both regimes t<-|t-statistic| and t>|t-statistic|.  The author of the post weaves in an explanation of 2 standard errors, but I suspect that this is unrelated to seeking the p-value for a specific t-statistic, especially since the value of |t| for 5% significance can depart nontrivially from 2.0 for small sample size.  Is this correct?
The above has also been posted to usenet groups matlab and stat.
The data to explore this problem is linked to from this tutorial.  The data file is located eggs.txt.
Here is the code that I used to explore the model:
clear variables

fnameDiary = [mfilename '.out.txt'];
if logical(exist(fnameDiary,'file'))
   diary off
   delete(fnameDiary)
end % if
diary(fnameDiary) % Also turns on diary

if0mean=true % Subtracts out the mean
ifminBias=false % std err normalizes by df rather than n,
                % where n=number of observations and df is n less the
                % number of parameters estimated (from vgxdisp)
nMaxLag=3 %#ok<NOPTS>
clf

tbChicEgg = readtable( 'eggs.txt', 'Delimiter', '\t' );

seriesNames={'chic','egg'};
varChicEgg = vgxset( 'Series', seriesNames, 'n',2 );
estVarChicEgg(nMaxLag).var = vgxset;

chicEgg = table2array(tbChicEgg(:,seriesNames));
if if0mean
   chicEgg = bsxfun( @minus, chicEgg, mean(chicEgg) );
end
chicEgg0 = chicEgg(1:nMaxLag,:);
chicEgg = chicEgg(1+nMaxLag:end,:);
yrs = table2array(tbChicEgg(1+nMaxLag:end,'year'));

subplot(1+nMaxLag,1,1);
plotyy( yrs,chicEgg(:,1) , yrs,chicEgg(:,2) );

for maxLag = 1:nMaxLag

   [est.spec est.stdErr est.LLF est.W] = vgxvarx( ...
      vgxset( varChicEgg, 'nAR',maxLag ), ...
      chicEgg, NaN, chicEgg0, ...
      'StdErrType', 'all' ...
   );

   fprintf('-------------------------\nmaxLag=%g\n',maxLag);
   estVarChicEgg(maxLag).var = est %#ok<NOPTS>

   subplot(1+nMaxLag,1,1+maxLag)
   plotyy(yrs,est.W(:,1),yrs,est.W(:,2))

   vgxdisp(est.spec,est.stdErr,'DoFAdj',ifminBias);

end

%--------------------------------------------------
diary off
% error('Stopping before return.');
return

 A: First of all, that answer on Matlab Central is really bad. It doesn't really answer the question, and it gives some bad advice. The formula they give computes the $p$-value for a two-sided test given a $t$-statistic and its degrees of freedom. The information they failed to give is how to obtain the error degrees of freedom dfe from the output of vgxvarx. And saying, don't bother about $p$-values just look at whether $|t| > 2$, ignores a host of problems, among other things how to deal with multiple tests. I'm afraid I cannot clear up all of it, but here's a partial answer.
1) You are right that in estimating the model, error degrees of freedom are lost, and in order to properly compute the $t$-statistic and to properly interpret the $t$-statistic, we would need to know how many error degrees of freedom that are. It appears that using the syntax
vgxdisp(EstSpec, EstStdErrors, 'DoFAdj', true)

you can obtain properly computed $t$-statistics, but the output still does not specify the degrees of freedom, so you still cannot properly interpret the result, and in particular, you cannot compute the $p$-value. Without wading through the internals of vgxvarx and vgxdisp, I don't see a way to determine these degrees of freedom. My recommendation: Ask the MathWorks, and don't accept a crappy answer.
If you know you have a very large number of data points, you can approximately assume that dfe $\rightarrow \infty$, in which case the $t$-distribution becomes the standard normal distribution, and you can approximately compute pValue = 2 * normcdf(-abs(t)).
Regarding the factor you mention, sqrt(df/n), no, this is not for correcting dfs, but to convert from the standard deviation associated with a mean to the standard error of that mean (and regression coefficient estimates can be seen as weighted means). This formula is not appropriate here.
2) Yes, the factor 2 comes from the fact that the $p$-value for a two-sided test is being computed. You are right that this is not related to the "rule of thumb" mentioned there, and you are right that this rule of thumb can be severely off. The rule of thumb is however related to the normal distribution mentioned above: 2 * normcdf(-2) is about 0.0455.
A: Thanks to A. Donda for the answer to part 2.  I was able to investigate my conjecture for part 1 in my original post for this question.  Yes, maximum likelihood estimate of the standard error for the coefficient uses a normalization of $n$ (the number of observations) whereas minimum bias estimate uses a normalization of $n$ less the number of coefficients.  I counted the number of coefficients from the vgxdisp output for the case of (maximum lag)=3.  I then compared the standard errors of the estimated coefficients with the vgxvarx's 'DoFAdj' set to true and to false (thanks, A. Donda), which yields the minimum bias and maximum likelihood estimates, respectively.
My comparison of the two sets of standard errors consisted of multiplying the maximum likelihood standard errors by $n$ to cancel out the normalization by $n$, then dividing by $n$ less the number of parameters (the normalization for minimum bias).  Note that I had it backward in my original post.  
There is a huge caveat though.  I had to increase the number of parameters by exactly 1 in order for the numbers to match.  So there must be an extra estimated parameter which isn't listed in the vgxdisp output.  Perhaps it is the cross-covariances between the noise sources for the two time series in the model.  It only shows up once in the vgxdisp output, but the cross-covariances were considered to be the same parameter because they should be identical.  Maybe for regression purposes, they are not.
Here is the code that converts maximum likelihood estimates of standard error to minimum bias estimates, with the number of estimated parameters $nEstParam$ being 1 more than the number listed in vgxdisp's output:
n=74; nEstParam=3*4+3+1; x=sqrt(n/(n-nEstParam)) * ...
[
   0.138092
    18.6127
 0.00102525
   0.138188
   0.151676
    29.8176
  0.0011261
   0.221377
   0.131599
    19.4701
0.000977039
   0.144553
5.81405e+07
     349742
    3204.78
]; for i=1:length(x); disp(x(i)), end

P.S. I realize that the solver is reporting poor conditioning of the matrix.  I am guessing that it is due to the apparent unit root in the 2nd time series, as hinted by the 1st subplot.
NOTE: As much as this answer confirms what is done in estimating the parameters, it still doesn't answer the question of whether the maximum likelihood or minimum variance estimate is better.  Not just better in a general sense (which probably hs no answer), but which is better for estimating p-value.  I am not sure whether there is a basis that can be put forth to favour one over the other.  Thoughts?
New info: I feel it is important to share that the mystery extra degree of freedom is not due to counting all 4 coefficients in the 2x2 covariance matrix, as suggested by the 3+1 contribution to nEstParam above. This was tested by selecting a diagonal covariance matrix (CovarType='diagonal' for vgxvar). We would expect that replacing 3+1 by 2 in nEstParam would yield the right conversion from ML stderr to min bias. It doesn't, and it turns out that the right conversion factor results from replacing 3+1 by 2+1 in nEstParam above. I suspect that vgxdisp is using the wrong degrees of freedom to calculate stderr because the resulting nEstParam doesn't even match the number of active parameters from vgxcount. I could be wrong, but until the extra DoF can be explained, it's the only possible conclusion. Hence, if I had to extract stderr from vgxdisp, I would request a ML stderr and correct it using the DoF from vgxcount. This is not ideal because it is still guessing, but it's the most sensible guess so far.
