%Create and SA fictitious series (CPS)

clear;
clc;
close all;
addpath(fullfile(cd, '/iris'))
irisstartup

jstr={'U1.spc','U2.spc','U3.spc','U4.spc','E1.spc','E2.spc','E3.spc','E4.spc'};

m=xlsread('cps.xls','nsa','C2:J449');
for j=1:4;  
    x=tseries(mm(1976,1):mm(2013,4),m(:,j));
    xsa=x12(x,Inf,'specfile=',char(jstr(j)));
    xsa=double(xsa);
    msa(:,j)=xsa;
end;
smsa=sum(msa')';
for j=1:11;
    param0(j,1)=(smsa(390)+((j/12)*(smsa(402)-smsa(390))))/smsa(390+j);
end;
param0=double(param0);


param=param0;
v(1)=myfuncu(param,m,jstr,[1;2;3;4]);
disp(v(1));
for k=1:100;
for j=1:11; 
    paramu=param; paramu(j)=paramu(j)+0.00001;
    paramd=param; paramd(j)=paramd(j)-0.00001;
    d(j,1)=(myfuncu(paramu,m,jstr,[1;2;3;4])-myfuncu(paramd,m,jstr,[1;2;3;4]))/0.00002;
end;
for j=1:101;
mn(j,1)=myfuncu(param-((1e-8)*(j-1)*d),m,jstr,[1;2;3;4]);
end;
[m1,m2]=min(mn);
param=param-((1e-8)*m2*d);
v(k+1)=myfuncu(param,m,jstr,[1;2;3;4]);
disp(param);
disp(v(k+1));
if v(k+1)<5; break; end;
end; %close the k=1 to 100 loop

clear msa;
m=xlsread('cps.xls','nsa','C2:J449');
m_alt=m(:,1:4);
for j=1:11; 
    m_alt(390+j,:)=m_alt(390+j,:)*param(j); 
end;
for j=1:4;  
    x=tseries(mm(1976,1):mm(2013,4),m_alt(:,j));
    xsa=x12(x,Inf,'specfile=',char(jstr(j)));
    xsa=double(xsa);
    msa(:,j)=xsa;
end;

xlswrite('cps.xls',msa,'sa_fict','C2:F449');

%-------------------------------------------------------%

m=xlsread('cps.xls','nsa','C2:J449');
for j=1:4;  
    x=tseries(mm(1976,1):mm(2013,4),m(:,j+4));
    xsa=x12(x,Inf,'specfile=',char(jstr(j+4)));
    xsa=double(xsa);
    msa(:,j)=xsa;
end;

smsa=sum(msa')';
for j=1:11;
    param0(j,1)=(smsa(390)+((j/12)*(smsa(402)-smsa(390))))/smsa(390+j);
end;
param0=double(param0);


param=param0;
v(1)=myfuncu(param,m,jstr,[5;6;7;8]);
disp(v(1));
for k=1:100;
for j=1:11; 
    paramu=param; paramu(j)=paramu(j)+0.00001;
    paramd=param; paramd(j)=paramd(j)-0.00001;
    d(j,1)=(myfuncu(paramu,m,jstr,[5;6;7;8])-myfuncu(paramd,m,jstr,[5;6;7;8]))/0.00002;
end;
for j=1:101;
mn(j,1)=myfuncu(param-((1e-10)*(j-1)*d),m,jstr,[5;6;7;8]);
end;
[m1,m2]=min(mn);
param=param-((1e-10)*m2*d);
v(k+1)=myfuncu(param,m,jstr,[5;6;7;8]);
disp(v(k+1));
if v(k+1)<5; break; end;
end; %close the k=1 to 100 loop

clear msa;
m=xlsread('cps.xls','nsa','C2:J449');
m_alt=m(:,5:8);
for j=1:11; 
    m_alt(390+j,:)=m_alt(390+j,:)*param(j); 
end;
for j=1:4;  
    x=tseries(mm(1976,1):mm(2013,4),m_alt(:,j));
    xsa=x12(x,Inf,'specfile=',char(jstr(j+4)));
    xsa=double(xsa);
    msa(:,j)=xsa;
end;

xlswrite('cps.xls',msa,'sa_fict','G2:J449');



irisfinish

