function [pi_hat,Y,X,Y_initial,Y_tilde_fit,err, Y_tilde, X_tilde]=VAR_LP(vardata,p,c,st_f)

[Y,X,Y_initial]=SUR(vardata,p,c);

Y_tilde=Y; for t=1:size(Y,1); Y_tilde(t,:)=Y(t,:)/st_f(t); end
X_tilde=X; for t=1:size(Y,1); X_tilde(t,:)=X(t,:)/st_f(t); end

pi_hat=(X_tilde'*X_tilde)\X_tilde'*Y_tilde;
Y_tilde_fit=X_tilde*pi_hat; %Fitted value of Y
err=Y_tilde-Y_tilde_fit; %Residuals

end