In[]:=
Clear["Global`*"];​​(****CreateDataforTesting*******)​​nn2=2;​​nbv=7;​​rr=10;​​data0=Table[{x1,x2,x3,x4,x5,x6,x7,​​1.2x1-3.4x2+10x3+2.8x4+1.6x5-1.9x6+2.8x7+​​RandomReal[rr]},{x1,RandomReal[rr,nn2]},{x2,​​RandomReal[rr,nn2]},{x3,RandomReal[rr,nn2]},{x4,​​RandomReal[rr,nn2]},{x5,RandomReal[rr,nn2]},{x6,​​RandomReal[rr,nn2]},{x7,RandomReal[rr,nn2]}];​​data0=Flatten[data0,nbv-1];​​data=Drop[data0,None,{nbv+1}];​​resp=Drop[data0,None,{1,nbv}];​​resp=Flatten[resp];​​(******Fitamodellm1*******)
In[]:=
lm1=LinearModelFit[{data,resp}]​​s1=(lm1["FitResiduals"])^2;
(*******Fitthesquarederrorstoamodelwiththesamepredictors\​lm2****)
In[]:=
lm2=LinearModelFit[{data,s1}];
(******ComputetheBreusch–Paganteststatistic*****)
bp=(Variance[s1](Length[data]-1)-Total[lm2["FitResiduals"]^2])/​​2/(Total[s1]/Length[data])^2​​(******OR*******)
bp2=With[{sqResids=​​lm1["FitResiduals"]^2},(Variance[sqResids](Length[data]-1)-​​Total[lm2["FitResiduals"]^2])/2/(Total[sqResids]/Length[data])^2]
(*********Computethep-value****)
dr=Length[lm1["BestFitParameters"]]-1;​​1-CDF[ChiSquareDistribution[dr],bp]