gausset;
library cmlmt,pgraph;
#include cmlmt.sdf

let maxH=21;
let Tboots=500000;
                                                    
/* The data is loaded                                                       */
load z[273,1]=EURUSD.txt;
true_z=z[253:273,1]; 
zAR=z[249:273,1];
z=z[251:273,1];

struct PV x0_GMAR;
x0_GMAR = pvCreate;
x0_GMAR = pvPacki(x0_GMAR,0.0429|-0.0116|1.2656|-0.2989,"varphi",1);
x0_GMAR = pvPacki(x0_GMAR,0.0581|0.0104,"sigma",2);
x0_GMAR = pvPacki(x0_GMAR,0.627,"alpha",3);

struct PV x0_GMAR_e;
x0_GMAR_e = pvCreate;
x0_GMAR_e = pvPacki(x0_GMAR_e,0.0407|-0.0111|1.2618|-0.2964,"varphi",1);
x0_GMAR_e = pvPacki(x0_GMAR_e,0.0583|0.0106,"sigma",2);
x0_GMAR_e = pvPacki(x0_GMAR_e,0.5862,"alpha",3);

struct PV x0_WL;
x0_WL = pvCreate;
x0_WL = pvPacki(x0_WL,0.0101|1.2570|-0.2715,"phi1",1);
x0_WL = pvPacki(x0_WL,0.0057,"phi2",2);
x0_WL = pvPacki(x0_WL,0.0563,"s1",3);
x0_WL = pvPacki(x0_WL,0.0092,"s2",4);
x0_WL = pvPacki(x0_WL,0.0328|2.4019,"beta",5);

struct PV x0_AR;
x0_AR = pvCreate;
x0_AR = pvPacki(x0_AR,0.0097|1.2783|-0.4187|0.3093|-0.1867,"alpha",1);
x0_AR = pvPacki(x0_AR,0.0366,"sigma",2);

proc f_ARp(struct PV x,z_d);
local a,s,forecast;
a = pvUnpack(x,1);
s = pvUnpack(x,2);
forecast=a[1]+a[2]*z_d[1]+a[3]*z_d[2]+a[4]*z_d[3]+a[5]*z_d[4];
retp(forecast);
endp;

proc f_GMAR(struct PV x, z_d);
local a,s,alpha,AA1,Gamma1,Gamma2,Ga1,Ga2,muu1,muu2;
a = pvUnpack(x,1);
s = pvUnpack(x,2);
alpha = pvUnpack(x,3);
AA1=((a[3]~a[4])|(1~0)); 
Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0);
Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0);
Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]);
Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]);
muu1=ones(2,1)*a[1]/(1-a[3]-a[4]);
muu2=ones(2,1)*a[2]/(1-a[3]-a[4]);
local ey1,ey2,forecast;
ey1=alpha*det(Ga1)^-0.5*exp(-(z_d-muu1)'*invpd(Ga1)*(z_d-muu1)/2);
ey2=(1-alpha)*det(Ga2)^-0.5*exp(-(z_d-muu2)'*invpd(Ga2)*(z_d-muu2)/2);
forecast=ey1/(ey1+ey2)*a[1]+ey2/(ey1+ey2)*a[2]+a[3]*z_d[1]+a[4]*z_d[2];
retp(forecast);
endp;

proc f_WLp(struct PV x, z_d);
local phi1,phi2,gamma1,gamma2,s1,s2,beta;
phi1 = pvUnpack(x,1);
phi2 = pvUnpack(x,2);
s1 = pvUnpack(x,3); 
s2 = pvUnpack(x,4);
beta = pvUnpack(x,5);
local apu,p,ey1,ey2,forecast;
apu=exp((1~z_d[2])*beta);
p=apu./(1+apu);
forecast=p*phi1[1]+(1-p)*phi2[1]+phi1[2:3]'*z_d[1:2];
retp(forecast);
endp;

/* 1-STEP CONDITIONAL EXPECTATION */

f_G=zeros(rows(true_z),1); e_G=zeros(rows(true_z),2); 
f_WL=zeros(rows(true_z),1); e_WL=zeros(rows(true_z),2);
f_G_e=zeros(rows(true_z),1); e_G_e=zeros(rows(true_z),2);
f_AR=zeros(rows(true_z),1); e_AR=zeros(rows(true_z),2);
for i (1, rows(true_z),1);
  f_G[i]=f_GMAR(x0_GMAR,z[2+i-1]|z[2+i-2]);
  f_WL[i]=f_WLp(x0_WL,z[2+i-1]|z[2+i-2]);
  f_G_e[i]=f_GMAR(x0_GMAR_e,z[2+i-1]|z[2+i-2]);
  f_AR[i]=f_ARp(x0_AR,zAR[4+i-1]|zAR[4+i-2]|zAR[4+i-3]|zAR[4+i-4]);
  e_G[i,1]=f_G[i]-true_z[i]; e_G[i,2]=(e_G[i,1])^2;
  e_WL[i,1]=f_WL[i]-true_z[i]; e_WL[i,2]=(e_WL[i,1])^2;
  e_G_e[i,1]=f_G_e[i]-true_z[i]; e_G_e[i,2]=(e_G_e[i,1])^2;
  e_AR[i,1]=f_AR[i]-true_z[i]; e_AR[i,2]=(e_AR[i,1])^2;
endfor;
/* AR(4) 1,...,MAXH CONDITIONAL EXPECTATION FORECASTS*/

f_ARH=zeros(rows(true_z),maxH); e_ARH=zeros(rows(true_z),2*maxH);
for j (1,maxH,1);
  for k (j,rows(true_z),1);
    if j==1;
      data=zAR[4+k-1]|zAR[4+k-2]|zAR[4+k-3]|zAR[4+k-4];
    elseif j==2;
      data=f_ARH[k-1,j-1]|zAR[4+k-2]|zAR[4+k-3]|zAR[4+k-4];
    elseif j==3;
      data=f_ARH[k-1,j-1]|f_ARH[k-2,j-2]|zAR[4+k-3]|zAR[4+k-4];
    elseif j==4;
      data=f_ARH[k-1,j-1]|f_ARH[k-2,j-2]|f_ARH[k-3,j-3]|zAR[4+k-4];
    else;
      data=f_ARH[k-1,j-1]|f_ARH[k-2,j-2]|f_ARH[k-3,j-3]|f_ARH[k-4,j-4];
    endif;  
    f_ARH[k,j]=f_ARp(x0_AR,data);
  endfor;
endfor;
for i (1,rows(true_z),1);
   e_ARH[i,1:maxH]=f_ARH[i,.]-true_z[i]; 
   e_ARH[i,(maxH+1):2*maxH]=(e_ARH[i,1:maxH])^2;
endfor;

screen on; print "AR(4) ready"; screen off;

/* 1,...,MAXH -STEP MC FORECASTS */
f_G_MC=zeros(rows(true_z),maxH); e_G_MC=zeros(rows(true_z),2*maxH); 
f_WL_MC=zeros(rows(true_z),maxH); e_WL_MC=zeros(rows(true_z),2*maxH);
f_G_e_MC=zeros(rows(true_z),maxH); e_G_e_MC=zeros(rows(true_z),2*maxH);
for k (1,rows(true_z),1);
  lth_a=minc(maxH|(rows(true_z)-k+1));
  a_f=zeros(Tboots,lth_a);  
  for i (1,Tboots,1);
    a_f[i,.]=f_GMAR22_MC(x0_GMAR,(z[2+k-1]|z[2+k-2]),lth_a)';   
  endfor;
  av_f=meanc(a_f);
  for j (1,lth_a,1);
     f_G_MC[j+k-1,j]=av_f[j];  
  endfor;
  screen on; print "GMAR(2,2) Conditional ready" k; screen off;
  for i (1,Tboots,1);
    a_f[i,.]=f_LMAR22_MC(x0_WL,(z[2+k-1]|z[2+k-2]),lth_a)';  
  endfor;
  av_f=meanc(a_f);
  for j (1,lth_a,1);
     f_WL_MC[j+k-1,j]=av_f[j];  
  endfor; 
  screen on; print "LMAR(2,2) ready" k; screen off;
  for i (1,Tboots,1);
    a_f[i,.]=f_GMAR22_MC(x0_GMAR_e,(z[2+k-1]|z[2+k-2]),lth_a)';   
  endfor;
  av_f=meanc(a_f);
  for j (1,lth_a,1);
     f_G_e_MC[j+k-1,j]=av_f[j];  
  endfor;  
  screen on; print "GMAR(2,2) Exact ready" k; print ""; screen off;
endfor;
for i (1,rows(true_z),1);
   e_G_MC[i,1:maxH]=f_G_MC[i,.]-true_z[i]; 
   e_G_MC[i,(maxH+1):2*maxH]=(e_G_MC[i,1:maxH])^2;
   e_WL_MC[i,1:maxH]=f_WL_MC[i,.]-true_z[i]; 
   e_WL_MC[i,(maxH+1):2*maxH]=(e_WL_MC[i,1:maxH])^2;
   e_G_e_MC[i,1:maxH]=f_G_e_MC[i,.]-true_z[i]; 
   e_G_e_MC[i,(maxH+1):2*maxH]=(e_G_e_MC[i,1:maxH])^2;
endfor;

proc f_GMAR22_MC(struct PV x0,z_d,T_out); 
local a,s,p,AA1,Gamma1,Gamma2,Ga1,Ga2,muu1,muu2;
a = pvUnpack(x0,1);
s = pvUnpack(x0,2);
p = pvUnpack(x0,3);
AA1=((a[3]~a[4])|(1~0)); 
Gamma1=inv(eye(2*2)-AA1.*.AA1)*(s[1]|0|0|0);
Gamma2=inv(eye(2*2)-AA1.*.AA1)*(s[2]|0|0|0);
Ga1=Gamma1[1:2]~(Gamma1[2]|Gamma1[1]);
Ga2=Gamma2[1:2]~(Gamma2[2]|Gamma2[1]);
muu1=ones(2,1)*a[1]/(1-a[3]-a[4]);
muu2=ones(2,1)*a[2]/(1-a[3]-a[4]);
local epsilon, boolean, data,ey1,ey2,bool;
epsilon=rndn(2+T_out,1);
data=zeros(2+T_out,1);
bool=zeros(2+T_out,1);
ey1=p*det(Ga1)^(-0.5)*exp(-(z_d-muu1)'*invpd(Ga1)*(z_d-muu1)/2);
ey2=(1-p)*det(Ga2)^(-0.5)*exp(-(z_d-muu2)'*invpd(Ga2)*(z_d-muu2)/2);
boolean=RNDBernoulli(ey1/(ey1+ey2));
bool[2,1]=boolean;
data[1,1]=z_d[2]; data[2,1]=z_d[1];
for k (3,T_out+2,1); 
  if boolean==1;
    data[k]=a[1]+a[3]*data[k-1]+a[4]*data[k-2]+sqrt(s[1])*epsilon[k];
  else;
    data[k]=a[2]+a[3]*data[k-1]+a[4]*data[k-2]+sqrt(s[2])*epsilon[k];
  endif;
  ey1=p*det(Ga1)^(-0.5)*exp(-(data[(k-1):k]-muu1)'*invpd(Ga1)*(data[(k-1):k]-muu1)/2);
  ey2=(1-p)*det(Ga2)^(-0.5)*exp(-(data[(k-1):k]-muu2)'*invpd(Ga2)*(data[(k-1):k]-muu2)/2);
  boolean=RNDBernoulli(ey1/(ey1+ey2));
  bool[k,1]=boolean;
endfor;
bool=bool[2:(rows(bool)-1),1];
retp(data[3:T_out+2]);
endp;

proc f_LMAR22_MC(struct PV x0,z_d,T_out); 
local phi1,phi2,gamma1,gamma2,s1,s2,beta;
phi1 = pvUnpack(x0,1);
phi2 = pvUnpack(x0,2);
s1 = pvUnpack(x0,3); 
s2 = pvUnpack(x0,4);
beta = pvUnpack(x0,5);
local epsilon, boolean, data, ey1,ey2,bool,apu,p;
epsilon=rndn(2+T_out,1);
data=zeros(2+T_out,1);
bool=zeros(2+T_out,1);
data[1:2]=z_d[2]|z_d[1];
apu=exp((1~data[1])*beta);
p=apu./(1+apu);
boolean=RNDBernoulli(p);
bool[2,1]=boolean;
for k (3,T_out+2,1); 
  if boolean==1;
    data[k]=phi1[1]+phi1[2:3]'*data[k-1:k-2]+sqrt(s1).*epsilon[k];
  else;
    data[k]=phi2[1]+phi1[2:3]'*data[k-1:k-2]+sqrt(s2).*epsilon[k];
  endif;
  apu=exp((1~data[k-1])*beta);
  p=apu./(1+apu);
  boolean=RNDBernoulli(p);
  bool[k,1]=boolean;
endfor;
bool=bool[2:(rows(bool)-1),1];
retp(data[3:T_out+2]);
endp;

proc RNDBernoulli(p);
local bool;
bool=rndu(1,1);
if p<0;   retp (error(0));
elseif p>1;   retp (error(0));
else;
  if bool<=p;    retp(1); 
  else;      retp(2);
  endif;
endif;
endp;
    
let pr=1; 
outwidth 120;
output file = forecasting_GMAR22_restricted_WL_rep500000.asc on;
screen off; 

print "";
print "GMAR conditional model parameters";
print pvGetParNames(x0_GMAR)'; 
print pvGetParVector(x0_GMAR)';
print "";
print "GMAR exact model parameters";
print pvGetParNames(x0_GMAR_e)'; 
print pvGetParVector(x0_GMAR_e)';
print "";
print "Wong and Li model parameters";
print pvGetParNames(x0_WL)'; 
print pvGetParVector(x0_WL)';
print "";
print "AR model parameters";
print pvGetParNames(x0_AR)'; 
print pvGetParVector(x0_AR)';
print "";
print "**************************************************************";
print "ONE-STEP FORECASTS ";
print "       True value       f_GMAR           f_GMAR_e         f_WL          f_AR";
for i (1, rows(true_z),1);
  print true_z[i]~f_G[i]~f_G_e[i]~f_WL[i]~f_AR[i];
endfor;
print "        (f_GMAR-tv)      (f_GMAR_e-tv)    (f_WL-tv)       (f_AR-tv)";
for i (1, rows(true_z),1);
  print e_G[i,1]~e_G_e[i,1]~e_WL[i,1]~e_AR[i,1];
endfor;
print "       (f_GMAR-tv)^2    (f_GMAR_e-tv)^2  (f_WL-tv)^2   (f_AR-tv)^2";
for  i (1, rows(true_z),1);
  print e_G[i,2]~e_G_e[i,2]~e_WL[i,2]~e_AR[i,2];
endfor;
print "";
GMAR_val=meanc(abs(e_G[.,1]))~meanc(e_G[.,2])~median(abs(e_G[.,1]))~median(e_G[.,2]);
print "GMAR model conditional";
print "      MAPE            MSPE             MedAPE           MedSPE"; 
print 100*GMAR_val./GMAR_val;
print "";
print 100*(meanc(abs(e_G_e[.,1]))~meanc(e_G_e[.,2])~median(abs(e_G_e[.,1]))~median(e_G_e[.,2]))./GMAR_val;
print "";
print 100*(meanc(abs(e_WL[.,1]))~meanc(e_WL[.,2])~median(abs(e_WL[.,1]))~median(e_WL[.,2]))./GMAR_val;
print "";
print 100*(meanc(abs(e_AR[.,1]))~meanc(e_AR[.,2])~median(abs(e_AR[.,1]))~median(e_AR[.,2]))./GMAR_val;
print "**************************************************************";
print "*************";
print "MC FORECASTS ";
print "*************";
print "**************************************************************";
for k (1,maxH,1);
  print "";
  print "step length in forecast: " k;
  print "       True value       f_GMAR           f_GMAR_e         f_WL          f_AR_H";
  for i (k,rows(true_z),1);
    print true_z[i]~f_G_MC[i,k]~f_G_e_MC[i,k]~f_WL_MC[i,k]~f_ARH[i,k];
  endfor;
  print "        (f_GMAR-tv)      (f_GMAR_e-tv)    (f_WL-tv)       (f_AR_H-tv)";
  for i (k,rows(true_z),1);
    print e_G_MC[i,k]~e_G_e_MC[i,k]~e_WL_MC[i,k]~e_ARH[i,k];
  endfor;
  print "       (f_GMAR-tv)^2    (f_GMAR_e-tv)^2  (f_WL-tv)^2   (f_AR_H-tv)^2";
  for i (k,rows(true_z),1);
    print e_G_MC[i,maxH+k]~e_G_e_MC[i,maxH+k]~e_WL_MC[i,maxH+k]~e_ARH[i,maxH+k];
  endfor;
endfor;

MAPEs_GMARc=zeros(maxH,4); MAPEs_GMARe=zeros(maxH,4); MAPEs_LMAR=zeros(maxH,4); MAPEs_AR=zeros(maxH,4);
for k (1,maxH,1);  
  print "";
  GMAR_val=meanc(abs(e_G_MC[k:rows(true_z),k]))~meanc(e_G_MC[k:rows(true_z),maxH+k])~median(abs(e_G_MC[k:rows(true_z),k]))~median(e_G_MC[k:rows(true_z),maxH+k]);
  print "step length in forecast: " k;
  print "";
  print "      MAPE            MSPE             MedAPE           MedSPE"; 
  MAPEs_GMARc[k,.]=100*GMAR_val./GMAR_val;
  print MAPEs_GMARc[k,.];
  print "";
  MAPEs_GMARe[k,.]=100*(meanc(abs(e_G_e_MC[k:rows(true_z),k]))~meanc(e_G_e_MC[k:rows(true_z),maxH+k])~median(abs(e_G_e_MC[k:rows(true_z),k]))~median(e_G_e_MC[k:rows(true_z),maxH+k]))./GMAR_val;
  print MAPEs_GMARe[k,.];
  print "";
  MAPEs_LMAR[k,.]=100*(meanc(abs(e_WL_MC[k:rows(true_z),k]))~meanc(e_WL_MC[k:rows(true_z),maxH+k])~median(abs(e_WL_MC[k:rows(true_z),k]))~median(e_WL_MC[k:rows(true_z),maxH+k]))./GMAR_val;
  print MAPEs_LMAR[k,.];
  print "";
  MAPEs_AR[k,.]=100*(meanc(abs(e_ARH[k:rows(true_z),k]))~meanc(e_ARH[k:rows(true_z),maxH+k])~median(abs(e_ARH[k:rows(true_z),k]))~median(e_ARH[k:rows(true_z),maxH+k]))./GMAR_val;
  print MAPEs_AR[k,.];
  print "**************************************************************";
endfor;  

output off;

outwidth 120;
output file = forecasting_errors_rep500000.asc on;
screen off; 
print "GMAR conditional";
print MAPEs_GMARc;
print "GMAR exact";
print MAPEs_GMARe;
print "LMAR";
print MAPEs_LMAR;
print "AR ";
print MAPEs_AR;

output off;
screen on;

apuria=seqa(1,1,12);
graphset;
 ylabel("");          
 xlabel("Step length"); 
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4)";
 _ptek="Forecast_MAPE12.tkf";
 title("MAPE of forecasts");
 xy(apuria,MAPEs_GMARc[1:12,1]~MAPEs_GMARe[1:12,1]~MAPEs_LMAR[1:12,1]~MAPEs_AR[1:12,1]);
 ret = tkf2eps("Forecast_MAPE12.tkf","Forecast_MAPE12.eps");

graphset;
 ylabel("");          
 xlabel("Step length"); 
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4)";
 _ptek="Forecast_MSPE12.tkf";
 title("MSPE of forecasts");
 xy(apuria,MAPEs_GMARc[1:12,2]~MAPEs_GMARe[1:12,2]~MAPEs_LMAR[1:12,2]~MAPEs_AR[1:12,2]);
 ret = tkf2eps("Forecast_MSPE12.tkf","Forecast_MSPE12.eps");

graphset;
 ylabel("");          
 xlabel("Step length"); 
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4)";
 _ptek="Forecast_MedAPE12.tkf";
 title("MedAPE of forecasts");
 xy(apuria,MAPEs_GMARc[1:12,3]~MAPEs_GMARe[1:12,3]~MAPEs_LMAR[1:12,3]~MAPEs_AR[1:12,3]);
 ret = tkf2eps("Forecast_MedAPE12.tkf","Forecast_MedAPE12.eps");


graphset;
 ylabel("");          
 xlabel("Step length");
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4)";
 xtics(0,12,4,2);
 ytics(90,230,40,20);
 _ptek="Forecast_MedSPE12.tkf";
 title("MedSPE of forecasts");
 xy(apuria,MAPEs_GMARc[1:12,4]~MAPEs_GMARe[1:12,4]~MAPEs_LMAR[1:12,4]~MAPEs_AR[1:12,4]);
 ret = tkf2eps("Forecast_MedSPE12.tkf","Forecast_MedSPE12.eps");
 
apuria=seqa(2010+1/13,1/12,(rows(true_z)));
forecasts=zeros(1,4);
for k (1,rows(true_z),1);
  forecasts=forecasts|(f_G_MC[k,k]~f_G_e_MC[k,k]~f_WL_MC[k,k]~f_ARH[k,k]);
endfor;
forecasts=forecasts[2:rows(forecasts),.];
graphset;
 ylabel("");          
 xlabel("Time"); 
 _pdate="";            
 _pltype ={1,2,3,4,6};
 _pcolor={0,0,0,0,0};
 _plwidth={1,2,1,2,1};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4) \0"\    
          "True value";
 _ptek="steps_1_21.tkf";
 title("Forecasts: steps 1 to 21");
 xy(apuria,forecasts~true_z);
 ret = tkf2eps("steps_1_21.tkf","steps_1_21.eps");

for k (1,minc(maxH|20),1);  
 nimi_t=ftos(k,"%*.*lf step MC predictions",1,0);
 nimitkf=ftos(k,"Forecasts_step_%*.*lf.tkf",1,0);
 nimieps=ftos(k,"Forecasts_step_%*.*lf.eps",1,0);
 apuria=seqa(2010+k/13,1/12,(rows(true_z)-k+1)); 
 graphset;
 ylabel("");          
 xlabel("Time"); 
 _pdate="";            
 _pltype ={1,2,3,4,6};
 _pcolor={0,0,0,0,0};
 _plwidth={1,2,1,2,1};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4) \0"\    
          "True value";
 _ptek=nimitkf;
 title(nimi_t);
 xy(apuria,f_G_MC[k:rows(f_G_MC),k]~f_G_e_MC[k:rows(f_G_MC),k]~f_WL_MC[k:rows(f_G_MC),k]~f_ARH[k:rows(f_G_MC),k]~true_z[k:rows(f_G_MC)]);
 ret = tkf2eps(nimitkf,nimieps);

 nimi_t=ftos(k,"%*.*lf step MC predictions; errors",1,0);
 nimitkf=ftos(k,"F_errors_%*.*lf.tkf",1,0);
 nimieps=ftos(k,"F_errors_%*.*lf.eps",1,0);
 graphset;
 ylabel("");          
 xlabel("Time");
 _pdate="";             
 _pltype ={1,3,5,6};
 _pcolor={0,0,0,0};
 _plegctl= { 2, 5, 1, 1.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
           "GMAR exact \0"\
           "Wong & Li \0"\
           "AR(4) ";    
 _ptek=nimitkf;
 title(nimi_t);     
 xy(apuria,e_G_MC[k:rows(e_G_MC),k]~e_G_e_MC[k:rows(e_G_MC),k]~e_WL_MC[k:rows(e_G_MC),k]~e_ARH[k:rows(e_ARH),k]);
 ret = tkf2eps(nimitkf,nimieps);

 nimi_t=ftos(k,"%*.*lf step MC predictions; squared errors",1,0);
 nimitkf=ftos(k,"F_S_errors_%*.*lf.tkf",1,0);
 nimieps=ftos(k,"F_S_errors_%*.*lf.eps",1,0);
 graphset;
 ylabel("");          
 xlabel("Time");
 _pdate="";             
 _pltype ={1,3,5,6};
 _pcolor={0,0,0,0};
 _plegctl= { 2, 5, 0.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
           "GMAR exact  \0"\
           "Wong & Li  \0"\
           "AR(4)";    
 _ptek=nimitkf;
 title(nimi_t);     
 xy(apuria,e_G_MC[k:rows(e_G_MC),maxH+k]~e_G_e_MC[k:rows(e_G_MC),maxH+k]~e_WL_MC[k:rows(e_G_MC),maxH+k]~e_ARH[k:rows(e_G_MC),maxH+k]);
 ret = tkf2eps(nimitkf,nimieps);
endfor;

MAPEs_GMARc=100*ones(21,4);

apuria=seqa(1,1,12);
graphset;
 ylabel("");          
 xlabel("Step length"); 
 xtics(0,12,2,2);
 ytics(95,115,5,0);
 _pdate="";   
 _pstype ={1};         
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 5, 1.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "Wong & Li \0"\ 
          "AR(4)";
 _ptek="Forecast_MAPE12.tkf";
 title("MAPE of forecasts");
 xy(apuria,MAPEs_GMARc[1:12,1]~MAPEs_GMARe[1:12,1]~MAPEs_LMAR[1:12,1]~MAPEs_AR[1:12,1]);
 ret = tkf2eps("Forecast_MAPE12.tkf","Forecast_MAPE12.eps");

graphset;
 ylabel("");          
 xlabel("Forecast horizon"); 
 xtics(0,12,2,2);
 ytics(90,130,10,0); 
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 5, 1.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "LMAR \0"\ 
          "AR(4)";
 _ptek="Forecast_MSPE12.tkf";
 title("");
 xy(apuria,MAPEs_GMARc[1:12,2]~MAPEs_GMARe[1:12,2]~MAPEs_LMAR[1:12,2]~MAPEs_AR[1:12,2]);
 ret = tkf2eps("Forecast_MSPE12.tkf","Forecast_MSPE12.eps");

graphset;
 ylabel("");          
 xlabel("Forecast horizon");  
 xtics(0,12,2,2);
 ytics(90,150,10,0);
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 2.75, 4 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "LMAR \0"\ 
          "AR(4)";
 _ptek="Forecast_MedAPE12.tkf";
 title("");
 xy(apuria,MAPEs_GMARc[1:12,3]~MAPEs_GMARe[1:12,3]~MAPEs_LMAR[1:12,3]~MAPEs_AR[1:12,3]);
 ret = tkf2eps("Forecast_MedAPE12.tkf","Forecast_MedAPE12.eps");


graphset;
 ylabel("");          
 xlabel("Forecast horizon"); 
 xtics(0,12,2,2);
 ytics(90,230,20,0);
 _pdate="";            
 _pltype ={6,1,3,6};
 _pcolor={8,0,0,0};
 _plwidth={4,4,4,4};
 _plegctl= { 2, 5, 2.75, 3.5 };  /* legend size and locations */
 _plegstr= "GMAR conditional \0"\
          "GMAR exact \0"\
          "LMAR \0"\ 
          "AR(4)";
 _ptek="Forecast_MedSPE12.tkf";
 title("");
 xy(apuria,MAPEs_GMARc[1:12,4]~MAPEs_GMARe[1:12,4]~MAPEs_LMAR[1:12,4]~MAPEs_AR[1:12,4]);
 ret = tkf2eps("Forecast_MedSPE12.tkf","Forecast_MedSPE12.eps");