library optmum, pgraph; optset; graphset;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Loading the data set @

load data[38,2] = c:\Eco7425\t22p7.txt;
c=data[.,1];
y=data[.,2];
n=rows(data);
r=4;

@ Initializing global variables @

beta=zeros(2,1); theta=zeros(2,1);

@ xi - Parameters over which NLS is performed
  Note optimization will be done by choosing values for xi and not beta and
  lambda @
@ The optimization command @

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@nonlinear least square (AR(2) model) to estimate beta1, beta2, theta1 and theta2@

xi=zeros(4,1);
grad=zeros(4,1);
stderr_xi=zeros(4,1);
{xi,ssqS,grad,retcode} = optmum(&ssqfunc,start());
sse_unrestricted=ssqS;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Calculating standard errors @

stderr_xi=se(xi);

@ Writing the output @
@ Recovering our parameters of interest by inverse transformation @

beta[1]=xi[1];
beta[2]=xi[2];
theta[1] =xi[3];
theta[2] =xi[4];
output file = c:\Eco7425\output\assignment2.out reset; output on;
"(b) Use the 38 observation and nonlinear least square (AR(2) model) to estimate beta1, beta2, theta1 and theta2";
"";
if (retcode==0);
"Computation Status: Normal convergence is achieved";
else;
"Computation Status: Abnormal termination";
endif;
"";
"Minimized sum of square of residual for AR(2) model =";; sse_unrestricted;
"";
"beta = "; beta;
"theta = "; theta;
"";
"Standard error of beta[1] = ";stderr_xi[1];
"Standard error of beta[2] = ";stderr_xi[2];
"Standard error of theta[1] = ";stderr_xi[3];
"Standard error of theta[2] = ";stderr_xi[4];
@"The gradient at the minimum point";grad;@
@"The return code =";retcode;@
output off;


@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@Testing marginal propensity to consume beta2=0.85@

xi=zeros(3,1);
grad=zeros(3,1);
stderr_xi=zeros(3,1);
{xi,ssqS,grad,retcode} = optmum(&ssqfunc1,start1());
sse_restricted=ssqS;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Calculating standard errors @

stderr_xi=se_1(xi);

j=1; k=n-2-r;
fx=((sse_restricted-sse_unrestricted)/j)/(sse_unrestricted/k);
fp=cdffc(fx,j,k);
output file = c:\Eco7425\output\assignment2.out; output on;
"";
"(c) Test the hypothesis that marginal propensity to consume equals 0.85(i.e., beta2=0.85) against the hypothesis that it does not.";
"";
if (retcode==0);
"Computation Status: Normal convergence is achieved";
else;
"Computation Status: Abnormal termination";
endif;
"";
"Minimized sum of square of residuals for model AR(2)(beta2=0.85) = ";; sse_restricted;
"";
"beta = "; beta;
"theta = "; theta;
"";
"";
"Standard error of beta[1] = ";stderr_xi[1];
"Standard error of theta[1] = ";stderr_xi[2];
"Standard error of theta[2] = ";stderr_xi[3];
"";
"The probability of Type I error, fp = ";fp;
"with the significance level of alpha = 0.05";
if (fp<=0.05);
"We can reject the hypothesis that marginal propensity to consume equals 0.85";
else;
"We can not reject the hypothesis that marginal propensity to consume equals 0.85";
endif;
output off;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@Testing whether AR(1) model would have been adequate, i.e., theta2=0@
xi=zeros(3,1);
grad=zeros(3,1);
stderr_xi=zeros(3,1);

{xi,ssqS,grad,retcode} = optmum(&ssqfunc2,start2());

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Calculating standard errors @
stderr_xi=se_2(xi);
sse_ar1=ssqS;
j=1; k=n-2-r;
fx=((sse_ar1-sse_unrestricted)/j)/(sse_unrestricted/k);
fp=cdffc(fx,j,k);

output file = c:\Eco7425\output\assignment2.out; output on;
"";
"(d) Testing whether AR(1) model would have been adequate (i.e., theta2=0)";
"";
if (retcode==0);
"Computation Status: Normal convergence is achieved";
else;
"Computation Status: Abnormal termination";
endif;
"";
"Minimized sum of squares of residual for AR(1) model = ";; sse_ar1;
"";
"Standard error of beta[1] = ";stderr_xi[1];
"Standard error of beta[2] = ";stderr_xi[2];
"Standard error of theta[1] = ";stderr_xi[3];
"";
"The probability of Type I error, fp = ";fp;
"with the significance level of alpha = 0.05";
if (fp<=0.05);
"We can reject the hypothesis that AR(1) model would have been adequate";
else;
"We can not reject the hypothesis that AR(1) model would have been adequate";
endif;
output off;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@Using linear least square to estimate the beta1 and beta2@

_olsres=1;  @ instructs GAUSS to compute and store OLS residuals @

@ The GAUSS OLS command @
nr=n-2;
yr = zeros(nr,1);
xr = zeros(nr,1);
yr = c[3:n,1];
xr = y[3:n,1];
{a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11}=ols(0,yr,xr); "";
sse_linear=0;
i=1;
do while i<=nr;
sse_linear = sse_linear+(a10[i])^2;
i=i+1; 
endo;
output file = c:\Eco7425\output\assignment2.out; output on;
"";
"(e) Use last 36 observation and linear least square to estimate beta1, beta2";
"";
"The minumu sum of square of residuals for linear model = ";sse_linear;
"";
"beta = "; a3;
"";
"Standard error of beta[1] = ";a6[1];
"Standard error of beta[2] = ";a6[2];
output off;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@Testing hypothesis theta1=0 and theta2=0@
j=2; k=n-2-r;
fx=((sse_linear-sse_unrestricted)/j)/(sse_unrestricted/k);
fp=cdffc(fx,j,k);
output file = c:\Eco7425\output\assignment2.out; output on;
"";
"(f) Testing hypothesis theta1=0 and theta2=0";
"";
"Minimized sum of squares of residuals for linear model = ";; sse_linear;
"";
"The probability of Type I error, fp = ";fp;
"with the significance level of alpha = 0.05";
if (fp<=0.05);
"We can reject the hypothesis that theta1=0 and theta2=0";
else;
"We can not reject the hypothesis that theta1=0 and theta2=0";
endif;
output off;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@Use linear least square to estimate the pis in the equation@

_olsres=1;  @ instructs GAUSS to compute and store OLS residuals @
nr=n-2;
yr = zeros(nr,1);
xr = zeros(nr,5);
yr = c[3:n,1];
xr[.,1]= c[2:n-1,1];
xr[.,2]= c[1:n-2,1];
xr[.,3]= y[3:n,1];
xr[.,4]= y[2:n-1,1];
xr[.,5]= y[1:n-2,1];

{a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11}=ols(0,yr,xr); "";
sse_linear_pi=0;
i=1;
do while i<=nr;
sse_linear_pi = sse_linear_pi+(a10[i])^2;
i=i+1; 
endo;

j=2; k=n-2-(r+2);
fx=((sse_unrestricted-sse_linear_pi)/j)/(sse_linear_pi/k);
fp=cdffc(fx,j,k);

output file = c:\Eco7425\output\assignment2.out; output on;
"";
"(g) Using linear least square to estimate the pis in the equation and"; 
"test the hypothesis that linear equation in (g) is"; 
"equavalent to the equation in (a) with certain restrictions";
"";
"The minumu sum of square of residuals for linear model equivalent to AR(2) = ";sse_linear_pi;
"";
"pi = "; a3;
"";
"Standard error of pi[1] = ";a6[1];
"Standard error of pi[2] = ";a6[2];
"Standard error of pi[3] = ";a6[3];
"Standard error of pi[4] = ";a6[4];
"Standard error of pi[5] = ";a6[5];
"Standard error of pi[6] = ";a6[6];
"";
"The probability of Type I error, fp = ";fp;
"with the significance level of alpha = 0.05";
if (fp<=0.05);
"We can reject the hypothesis that linear equation in (g) is equavalent to the equation in (a) with certain restrictions";
else;
"We can not reject the hypothesis that linear equation in (g) is equavalent to the equation in (a) with certain restrictions";
endif;
output off;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
proc start();
local ini_para;

beta[1]=0;
beta[2]=0.8;  
theta[1]=0.5;
theta[2] =0.5;

@ Theta values are chosen as the starting values for xi, the parameters
  over which optimization is performed @

ini_para=zeros(4,1);
ini_para[1]=beta[1];
ini_para[2]=beta[2];
ini_para[3]=theta[1];
ini_para[4]=theta[2];
retp(ini_para); 
endp;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Procedure to evaluate the sum of squares - the criterion function
  that is to be minimized @

proc ssqfunc(xi);
local eY,i,ssqs,temp;

beta[1]=xi[1];
beta[2]=xi[2];
theta[1]=xi[3];
theta[2]=xi[4];

@ Calculating the sum of squares @

eY=zeros(rows(data),1);
i=3; 
do while i<=rows(eY);
temp = beta[1]*(1-theta[1]-theta[2])+theta[1]*c[i-1]+theta[2]*c[i-2];
temp = temp + beta[2]*(y[i]-theta[1]*y[i-1]-theta[2]*y[i-2]);
eY[i] = (c[i]-temp)^2;
i=i+1; 
endo;
ssqS = sumc(eY);
retp(ssqs); 
endp;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Calculating standard errors for function ssqfunc@

proc se(xi);
local i,info,se,temp,sgmsq;
optset;
sgmsq=ssqfunc(xi)/rows(data);
info=gradcd(&se1,xi);
info=invswp(info); 
se=sqrt(2*sgmsq*diag(info));
retp(se); 
endp;

proc se1(xi);
local info;
info=gradcd(&ssqfunc,xi);
retp(info'); 
endp;

@!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!@
@The following procedure is for testing@

proc start1();
local ini_para;
beta[1]=0;
beta[2]=0.85;   
theta[1]=0.5;
theta[2] =0.5;
@ Theta values are chosen as the starting values for xi, the parameters
  over which optimization is performed @

ini_para=zeros(3,1);
ini_para[1]=beta[1];
ini_para[2]=theta[1];
ini_para[3]=theta[2];
retp(ini_para); 
endp;

proc ssqfunc1(xi);
local eY,i,ssqs,temp;
beta[1]=xi[1];
beta[2]=0.85;
theta[1]=xi[2];
theta[2]=xi[3];

@ Calculating the sum of squares @

eY=zeros(rows(data),1);
i=3; 
do while i<=rows(eY);
temp = beta[1]*(1-theta[1]-theta[2])+theta[1]*c[i-1]+theta[2]*c[i-2];
temp = temp + beta[2]*(y[i]-theta[1]*y[i-1]-theta[2]*y[i-2]);
eY[i] = (c[i]-temp)^2;
i=i+1; 
endo;

ssqS = sumc(eY);
retp(ssqs); 
endp;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Calculating standard errors for function ssqfunc1@

proc se_1(xi);
local i,info,se,temp,sgmsq;
optset;
sgmsq=ssqfunc1(xi)/rows(data);
info=gradcd(&se1_1,xi);
info=invswp(info); 
se=sqrt(2*sgmsq*diag(info));
retp(se); 
endp;

proc se1_1(xi);
local info;
info=gradcd(&ssqfunc1,xi);
retp(info'); 
endp;


proc start2();
local ini_para;
beta[1]=0;
beta[2]=0.8;   
theta[1]=0.5;
theta[2] =0;
@ Theta values are chosen as the starting values for xi, the parameters
  over which optimization is performed @

ini_para=zeros(3,1);
ini_para[1]=beta[1];
ini_para[2]=beta[2];
ini_para[3]=theta[1];
retp(ini_para); 
endp;

proc ssqfunc2(xi);
local eY,i,ssqs,temp;
beta[1]=xi[1];
beta[2]=xi[2];
theta[1]=xi[3];
theta[2]=0;

@ Calculating the sum of squares @

eY=zeros(rows(data),1);
i=3; 
do while i<=rows(eY);
temp = beta[1]*(1-theta[1]-theta[2])+theta[1]*c[i-1]+theta[2]*c[i-2];
temp = temp + beta[2]*(y[i]-theta[1]*y[i-1]-theta[2]*y[i-2]);
eY[i] = (c[i]-temp)^2;
i=i+1; 
endo;
ssqS = sumc(eY);
retp(ssqs); 
endp;

@!!!!!!!!!!!!!!!!!!!!!!!!!@
@ Calculating standard errors for function ssqfunc2@

proc se_2(xi);
local i,info,se,temp,sgmsq;
optset;
sgmsq=ssqfunc2(xi)/rows(data);
info=gradcd(&se1_2,xi);
info=invswp(info); 
se=sqrt(2*sgmsq*diag(info));
retp(se); 
endp;

proc se1_2(xi);
local info;
info=gradcd(&ssqfunc2,xi);
retp(info'); 
endp;


end;