% IV for Sample 2 in Brady (2011)

load Nor_Ca %my prepared data loaded from the Matlab workspace (previously saved as a .mat file)

T = 89;  %Time in months for each region
nc = 23; %The number of regions

reg = data(:,16); %County numbers
dt = data(:,2); %date
y = data(:,3); %housing prices (natural log)
y1 = data(:,46); %housing prices lag
x1 = data(:,4);  %unemployment
x2 = data(:,13);  %construction (natural log)
x3 = data(:,43);  %realrate
x4 = data(:,44);  %population (natural log)
x5 = data(:,45); %industrial production (natural log)
x6 = data(:,17:19);  %season dummies
x7 = data(:,21:42);  %county dummies (22 for 23 counties)

%lags of the regressors
y2 = data(:,47); y3 = data(:,64); y4 = data(:,65); %more housing lags
x1_1 = data(:,50); x1_2 = data(:,51); x1_3 = data(:,68); x1_4 = data(:,69); %unemployment lags
x2_1 = data(:,48); x2_2 = data(:,49); x2_3 = data(:,66); x2_4 = data(:,67); %construction lags
x3_1 = data(:,52); x3_2 = data(:,53); x3_3 = data(:,54); x3_4 = data(:,55); %realrate lags
x4_1 = data(:,60); x4_2 = data(:,61); x4_3 = data(:,62); x4_4 = data(:,63); %population lags
x5_1 = data(:,56); x5_2 = data(:,57); x5_3 = data(:,58); x5_4 = data(:,59); %IP lags
yleads = leads(:,4:39);  %Housing price leads +1 through +36 (the third column in leads is the current hprice)

% 2) Here I define the initial weighting matrix.  James Lesage does provide
% code for automating this part, but I chose to define the matrix by
% "hand."

W=[[1,1,0,0,0,1,0,0,0,0,0,0,0,1,1,1,1,0,0,0,1,0,0;1,1,0,0,0,1,0,0,1,0,0,1,0,1,1,0,0,0,1,1,0,0,0;0,0,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0;0,0,0,1,1,0,1,1,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0;0,0,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0;1,1,0,0,0,1,0,0,1,0,0,0,0,1,0,0,0,0,1,1,0,0,0;0,0,0,1,1,0,1,0,0,0,0,0,1,0,0,0,1,0,0,0,1,0,0;0,0,0,1,0,0,0,1,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0;0,1,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,1,1,0,0,1;0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0;0,0,1,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0;0,1,1,0,0,0,0,0,0,0,1,1,0,0,1,0,0,0,1,0,0,0,1;0,0,0,1,0,0,1,1,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0;1,1,0,0,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0;1,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1,0,0,0,1,0,0;1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,1,0,0,0,0,0;1,0,0,0,0,0,1,1,0,0,0,0,1,0,1,1,1,1,0,0,1,0,0;0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,1,1,1,0,0,0,0,0;0,1,0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,1,1,0,0,1;0,1,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0;1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0;0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0;0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,1;]];


x = zeros([23 23]);

for i=1:23
    for j=1:23
        if j~=i
    x(i,j)=W(i,j);
        end
    end
end



% In the next few lines I sort the data
lags = [y2 y3 y4 x1_1 x1_2 x1_3 x1_4 x2_1 x2_2 x2_3 x2_4 x3_1 x3_2 x3_3 x3_4 x4_1 x4_2 x4_3 x4_4 x5_1 x5_2 x5_3 x5_4];
all = [dt reg y y1 x1 x2 x3 x4 x5 x6 x7 yleads lags];
[n,m] = size(all);

% Here I rid the data matrix "all" of the observations lost due to the lags
% used in estimation (i.e., if you decide on an AR(1), or AR(2) . . . )
for k = 1:n;
     %if k == 2671, break, end  %This should be plus 1 of the final sample (after taking lags)
     if k == 2048
         break, end
 if (all(k,4) == 0.0012); %Since I sorted the data in Excel, I had some zeros (after taking lags), which I changed to .0012 to import into Matlab
    all(k,:) = [];        %This line gets rid of rows with any cell equal to 0.0012.
 end;
end;

all = sortrows(all,1); %This sorts the matrix according to data ("dt"), this is required by Lesage's code used below.
[n,m] = size(all); 

y =  (all(:,3))*100; %Here I just redefine the variables in the "all" matrix
y1 = (all(:,4))*100;
x1 = all(:,5);
x2 = (all(:,6))*100;
x3 = all(:,7);
x4 = (all(:,8))*100;
x5 = (all(:,9))*100;
x6 = all(:,10:12);
x7 = all(:,13:34);
dt = all(:,1);
yleads = (all(:,35:70))*100; %36 Leads
lags = (all(:,71:93))*100;
W = x; %Rename the weighting matrix consistent with Lesage's code.

% Lesage's code: "sparse" "normw" and "slag" are from Lesage's spatial toolbox (see his webiste).  
WW = sparse(W);
W1 = normw(WW);
W = slag(W1,2);

% This code creates the "spatial lag" variable denoted in the literature as Wy.
% I also estimate the "lags" for the regressors for use in IV.
% These commands are based on Paul Elhorst's "sar_panel.m" file
% provided in Lesage's toolkit.

for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wy([t1:t2],1)= sparse(W)*y([t1:t2],1);
end
%Instruments
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx1([t1:t2],1)= sparse(W)*x1([t1:t2],1);
end
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx2([t1:t2],1)= sparse(W)*x2([t1:t2],1);
end
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx4([t1:t2],1)= sparse(W)*x4([t1:t2],1);
end
for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx1_1([t1:t2],1)= sparse(W)*lags([t1:t2],4);
end

for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wx2_1([t1:t2],1)= sparse(W)*lags([t1:t2],8);
end

for t=1:T
    t1=1+(t-1)*nc;t2=t*nc;
    Wy1([t1:t2],1)= sparse(W)*y1([t1:t2],1);
end

% Note: There are plenty of OLS or TSLS matlab programs out there (see Lesage's
% collection of code, for example). I include my own code *in* the program
% for my own satisfaction (as a way to keep track of things).  
% You could easily cut lines 184 through 270 and go straigt to the loop
% below (which uses Lesage's hwhite(y,X) program).  

% 3) OLS: 

% define nt-by-nt weithting matrix

W_nt=kron(eye(89),W);


% OLS with spatial regressor Shulin Shen write the OLS code
[n,m] = size(y);
cc = ones(n,1);
X = [cc Wy y1 x1 x2 x3 x4 x5 x6 x7];
[n,k] = size(X);
bhat = inv(X'*X)*(X'*y);  
ybar = mean(y);
sst=(y-ybar)'*(y-ybar); % Shulin Shen add the total sum of squares
e_fs = y - X*bhat;  
y_OLS = X*bhat;
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(X'*X);  
se = sqrt(diag(varco));
tstat = bhat./se;
pvalue = 2*(1 - tcdf(abs(tstat),n-k));  
OLS_results.y_OLS = y_OLS;
OLS_results.bhat = bhat;
OLS_results.rho = bhat(3,:);
OLS_results.sse=sse;
OLS_results.e_fs = e_fs;
OLS_results.se = se;
OLS_results.tstat = tstat;
OLS_results.pvalue = pvalue;
OLS_results.r2=1-(n-k)*sse/sst; % Shulin Shen add the R-square
OLS_results.ad_r2=1-(n-1)/(n-k-2)*(n-k)*sse/sst; % Shulin Shen add the adjusted R-square

% Shulin Shen's White's Corrected Standard Errors

sigma=diag(e_fs.^2);

varco_white_beta=inv(X'*X)*(X'*sigma*X)*inv(X'*X);

se_white=sqrt(diag(varco_white_beta));

tstat_white=bhat./se_white;

OLS_results.pvalue_white=2*(1 - tcdf(abs(tstat_white),n-k));

OLS_RE=[OLS_results.bhat se_white OLS_results.pvalue_white OLS_results.pvalue OLS_results.tstat OLS_results.se];


%White's Corrected Standard Errors (from Lesage's Toolbox)
% white_se = hwhite(y,X);
% white_tstat = white_se.tstat;
%calclute the white se's
% se_white = bhat./white_tstat;
% pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));
% OLS_results.se_white = se_white;
% OLS_results.pvalue_white= pvalue_white;
% end of OLS 

%***************************************************************
%LM Test on residuals for the SAR model
%this program is pasted in from Lesage's program and modified for panel
%data

T=89; % number of time periods: 90 - one lag
N=23; % number of regions
nobs=N*T;

es = OLS_results.e_fs; %Residuals from the SAR-OLS estimation
rho = OLS_results.rho;
% recover variance of rho
rhot = OLS_results.tstat(2,1);
sige = OLS_results.sse;
stdt = rhot/rho;
stdt = 1/stdt;
varr = stdt*stdt;
[n junk] = size(W);
A = speye(N) - rho*sparse(W);
AI = inv(A);
W2 = sparse(W);
%W3=W2.*W2;
T22 =(T-1)* trace(W2*W2 + W2'*W2);
T21 = (T-1)*trace(W2'*W*AI + W2*W*AI);

lm1=0;
for t=1:T  %This "for" statement modifies Lesage's code for panel data; I used template from Elhorst's sar_panel.m
    t1=1+(t-1)*N;t2=t*N;
lm1 =lm1+ (es([t1:t2],:)'*W*es([t1:t2],:))/sige;
end
Tterm = (T22 - T21*T21*varr);
TI = inv(Tterm);
lmerr = lm1*lm1*TI;
prob = 1-chis_prb(lmerr,1);

OLS_LMSARresult.meth = 'lmsar';
OLS_LMSARresult.lm = lmerr;
OLS_LMSARresult.prob = prob;
OLS_LMSARresult.chi1   = 6.635;
OLS_LMSARresult.nobs = nobs;
OLS_LMSARresult.nvar = k;


% 4) Estimating the Linear Projections: This code is a simple loop that
%    estimates the leads of the dependent variable on the regressors.  


[p,q]=size(all);
cc = ones(p,1);
nobs=N*T;
% Execute the estimation for the leads of the dependent variable
% First drop the last nlost observations from both the ylead and x matrix
% Note: The loop below is for only one regressor. Either copy the code for
% any additional regressor of interest, or just change the row number for
% the "beta" and "se" lines.  

% x = [cc y1 WyI x1 x2 x3 x5 x4 x6 x7];%IV for spatial regressor
xx = [cc y1 Wy x1 x2 x3 x5 x4 x6 x7]; %Basic OLS

% Note thw code below is written assuming that line 122 has been executed
for t=1:30 %Length of the horizon
    nt = N*(T-t);
    ylead = yleads(1:nt,t);
    xlead = xx(1:nt,:);
    spat_lead(t,1) = hwhite(ylead,xlead); 
end
for t=1:30
    irf(t,1) = spat_lead(t,1).beta(3,:); %Collect the coefficient on the third term in x
    se2(t,1) = spat_lead(t,1).beta(3,:)/spat_lead(t,1).tstat(3,:);
end

irf_se_plus2 = irf+(se2*2);
irf_se_minus2 = irf-(se2*2);
irf_spatial_withbands = [irf_se_plus2 irf irf_se_minus2];

%Send irf_spatial_withbands to an Excel file if you wish:
%xlswrite('irf_spatial.xls', irf_spatial_withbands)

%The plot produced below is simple; modify as you see fit.
% Create figure
%figure1 = figure;
% Create axes
%axes('Parent',figure1);
%box('on');
%hold('all');
% Create plot
%plot(irf_spatial_withbands,'DisplayName','irf_spatial(1:25,1)');

xmarker=linspace(1,30,30);
ymarker=irf_spatial_withbands;

plot(xmarker,irf_se_plus2,'k',xmarker,irf_se_plus2, 'ks','LineWidth',1,'MarkerSize',5);
axis([1 30 -1 1]); 
%plot(irf_spatial_withbands,irf_spatial_withbands,'*');
hold on;
plot(xmarker,irf,'k','LineWidth',1);
hold on;
plot(xmarker,irf_se_minus2,'k',xmarker,irf_se_minus2, 'ks','LineWidth',1,'MarkerSize',5);
hold on;
plot(xmarker,linspace(0,0,30),'LineWidth',1.5);


% 4) IV 

W_nt=kron(eye(89),W);
W1_nt=kron(eye(89),W1);

% First Stage estimation
[n,m] = size(Wy);
cc = ones(n,1);
WX = [cc Wx1 Wx2 Wx4 Wy1 y1 x1 x2 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat = inv(WX'*WX)*(WX'*Wy);  
Wybar = mean(Wy);
e_fs = Wy - WX*bhat;  
Wy_Instr = WX*bhat;
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se = sqrt(diag(varco));
tstat = bhat./se;
pvalue = 2*(1 - tcdf(abs(tstat),n-k)); 

% First stage F-test

RSSR=(Wy-Wybar)'*(Wy-Wybar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_pvalue=1-fcdf(F,k-1,n-k);
FirstStage_results.F_test=F;
FirstStage_results.F_test_pvalue=F_pvalue;
FirstStage_results.Wy_Instr = Wy_Instr;
FirstStage_results.bhat = bhat;
FirstStage_results.e_fs = e_fs;
FirstStage_results.se = se;
FirstStage_results.tstat = tstat;
FirstStage_results.pvalue = pvalue;
WyI = Wy_Instr;

First_IV=[bhat pvalue tstat se];


%Second Stage --------------------------------------------
[n,m] = size(y);
cc = ones(n,1);
X = [cc WyI y1 x1 x2 x3 x5 x4 x6 x7];
[n,k] = size(X);
bhat = inv(X'*X)*(X'*y);  
ybar = mean(y);
e = y - X*bhat;  
y_pred = X*bhat;
sse = e'*e/(n-k);  
varco = sse*inv(X'*X);  
se = sqrt(diag(varco));
tstat = bhat./se;
pvalue = 2*(1 - tcdf(abs(tstat),n-k));  

%__________________________________________________________________________
%White's Corrected Standard Errors (from Lesage's Toolbox)
white_se = hwhite(y,X);
white_tstat = white_se.tstat;

% Much of this code may be redundant or useless to some; I include it to
% provide a nice "printout" in the command window.

results.k = k;    
results.n = n;
results.X = X;
results.y_pred = y_pred;
results.sse = sse;
results.bhat = bhat;
results.rho = bhat(3,:);
results.e = e;
results.se = se;
results.tstat = tstat;
results.pvalue = pvalue;
%calclute the white se's
se_white = bhat./white_tstat;
resluts.pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));


IV_RE=[results.bhat se_white resluts.pvalue_white results.pvalue results.tstat results.se]


%***************************************************************
%LM Test on residuals for the SAR model
%this program is pasted in from Lesage's program and modified for panel
%data
es = results.e; %Residuals from the SAR-IV estimation
rho = results.rho;
% recover variance of rho
rhot = results.tstat(2,1);
sige = results.sse;
stdt = rhot/rho;
stdt = 1/stdt;
varr = stdt*stdt;
[n junk] = size(W);
A = speye(N) - rho*sparse(W);
AI = inv(A);
W2 = sparse(W);

T22 =(T-1)* trace(W2*W2 + W2'*W2);
T21 = (T-1)*trace(W2'*W*AI + W2*W*AI);

lm1=0;
for t=1:T  %This "for" statement modifies Lesage's code for panel data; I used template from Elhorst's sar_panel.m
    t1=1+(t-1)*N;t2=t*N;
lm1 =lm1+ (es([t1:t2],:)'*W*es([t1:t2],:))/sige;
end
Tterm = (T22 - T21*T21*varr);
TI = inv(Tterm);
lmerr = lm1*lm1*TI;
prob = 1-chis_prb(lmerr,1);

IV_LMSARresult.meth = 'lmsar';
IV_LMSARresult.lm = lmerr;
IV_LMSARresult.prob = prob;
IV_LMSARresult.chi1   = 6.635;
IV_LMSARresult.nobs = nobs;
IV_LMSARresult.nvar = k;






% 5) IV 2

W_nt=kron(eye(89),W);
W1_nt=kron(eye(89),W1);

% First Stage estimation
[n,m] = size(Wy);
cc = ones(n,1);

x1_1=lags(:,4);
x2_1=lags(:,8);

WX = [cc Wx1 Wx2 Wx4 Wy1 x1_1 x2_1 y1 x3 x5 x4 x6 x7]; 
% WX = [cc Wx1 Wx2 Wx4 Wy1 lags y1 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat_wy = inv(WX'*WX)*(WX'*Wy); 
Wy_Instr = WX*bhat_wy;
Wybar = mean(Wy);
e_fs = Wy - WX*bhat_wy;  
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se_wy= sqrt(diag(varco));
tstat_wy= bhat_wy./se_wy;
pvalue_wy= 2*(1 - tcdf(abs(tstat_wy),n-k)); 
RSSR=(Wy-Wybar)'*(Wy-Wybar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F_wy=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_wy_pvalue=1-fcdf(F_wy,k-1,n-k);

First_IV_Wy=[bhat_wy pvalue_wy tstat_wy se_wy];

[n,m] = size(x1);
cc = ones(n,1);
WX = [cc Wx1 Wx2 Wx4 Wy1 Wx1_1 Wx2_1 y1 x3 x5 x4 x6 x7]; 
% WX = [cc Wx1 Wx2 Wx4 Wy1 x1_1 x2_1 y1 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat_x1=inv(WX'*WX)*(WX'*x1);
x1_Instr = WX*bhat_x1;
x1bar = mean(x1);
e_fs = x1 - WX*bhat_x1;  
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se_x1= sqrt(diag(varco));
tstat_x1= bhat_x1./se_x1;
pvalue_x1= 2*(1 - tcdf(abs(tstat_x1),n-k)); 
RSSR=(x1-x1bar)'*(x1-x1bar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F_x1=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_x1_pvalue=1-fcdf(F_x1,k-1,n-k);

First_IV_x1=[bhat_x1 pvalue_x1 tstat_x1 se_x1];

[n,m] = size(x2);
cc = ones(n,1);
WX = [cc Wx1 Wx2 Wx4 Wy1 Wx1_1 Wx2_1 y1 x3 x5 x4 x6 x7];
% WX = [cc Wx1 Wx2 Wx4 Wy1 x1_1 x2_1 y1 x3 x5 x4 x6 x7]; %Shulin Shen excluded lags from the author's code
[n,k] = size(WX);
bhat_x2=inv(WX'*WX)*(WX'*x2); 
x2_Instr = WX*bhat_x2;
x2bar = mean(x2);
e_fs = x2 - WX*bhat_x2;  
sse = e_fs'*e_fs/(n-k);  
varco = sse*inv(WX'*WX);  
se_x2= sqrt(diag(varco));
tstat_x2= bhat_x2./se_x2;
pvalue_x2= 2*(1 - tcdf(abs(tstat_x2),n-k)); 
RSSR=(x2-x2bar)'*(x2-x2bar); % Restricted sum of squared residuals
USSR=e_fs'*e_fs;
F_x2=((RSSR-USSR)/(k-1))/(USSR/(n-k)); 
F_x2_pvalue=1-fcdf(F_x2,k-1,n-k);
 
First_IV_x2=[bhat_x2 pvalue_x2 tstat_x2 se_x2];

%Second Stage --------------------------------------------
[n,m] = size(y);
cc = ones(n,1);
X = [cc Wy_Instr y1 x1_Instr x2_Instr x3 x5 x4 x6 x7];
[n,k] = size(X);
bhat = inv(X'*X)*(X'*y);  
ybar = mean(y);
e = y - X*bhat;  
y_pred = X*bhat;
sse = e'*e/(n-k);  
varco = sse*inv(X'*X);  
se = sqrt(diag(varco));
tstat = bhat./se;
pvalue = 2*(1 - tcdf(abs(tstat),n-k));  

%__________________________________________________________________________
%White's Corrected Standard Errors (from Lesage's Toolbox)
white_se = hwhite(y,X);
white_tstat = white_se.tstat;

% Much of this code may be redundant or useless to some; I include it to
% provide a nice "printout" in the command window.

results.k = k;    
results.n = n;
results.X = X;
results.y_pred = y_pred;
results.sse = sse;
results.bhat = bhat;
results.rho = bhat(3,:);
results.e = e;
results.se = se;
results.tstat = tstat;
results.pvalue = pvalue;
%calclute the white se's
se_white = bhat./white_tstat;
resluts.pvalue_white = 2*(1 - tcdf(abs(white_tstat),n-k));


IV2_RE=[results.bhat se_white resluts.pvalue_white results.pvalue results.tstat results.se];


%***************************************************************
%LM Test on residuals for the SAR model
%this program is pasted in from Lesage's program and modified for panel
%data
es = results.e; %Residuals from the SAR-IV estimation
rho = results.rho;
% recover variance of rho
rhot = results.tstat(2,1);
sige = results.sse;
stdt = rhot/rho;
stdt = 1/stdt;
varr = stdt*stdt;
[n junk] = size(W);
A = speye(N) - rho*sparse(W);
AI = inv(A);
W2 = sparse(W);

T22 =(T-1)* trace(W2*W2 + W2'*W2);
T21 = (T-1)*trace(W2'*W*AI + W2*W*AI);

lm1=0;
for t=1:T  %This "for" statement modifies Lesage's code for panel data; I used template from Elhorst's sar_panel.m
    t1=1+(t-1)*N;t2=t*N;
lm1 =lm1+ (es([t1:t2],:)'*W*es([t1:t2],:))/sige;
end
Tterm = (T22 - T21*T21*varr);
TI = inv(Tterm);
lmerr = lm1*lm1*TI;
prob = 1-chis_prb(lmerr,1);

IV2_LMSARresult.meth = 'lmsar';
IV2_LMSARresult.lm = lmerr;
IV2_LMSARresult.prob = prob;
IV2_LMSARresult.chi1   = 6.635;
IV2_LMSARresult.nobs = nobs;
IV2_LMSARresult.nvar = k;


