%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% FILE: ss.m
%%
%% Stefan Stuidu miettii marraskuussa 2020 mita tehda 
%% syyskuu 2021-elokuu 2022 
%%
%% Han on paattanyt vaihtoehdot
%%
%% a1 = Opiskelijavaihto Kiinaan
%% a2 = Keskittyminen opintoihin
%% a3 = Opintoja ja ruokalahetin toita iltaisin
%% a4 = Vahan opintoja ja paatyo isossa kv. yrityksessa
%%
%% Seuraavat skenaariot ovat Stefanin mielesta keskeisia
%%
%% s1 = Korona riehuu, yliopisto etaopetuksessa ja maailma lamassa
%% s2 = Korona riehuu, yliopisto etaopetuksessa ja maailman talous elpyy
%% s3 = Korona riehuu, yliopisto normiopetuksessa ja maailma lamassa
%% s4 = Korona riehuu, yliopisto normiopetuksessa ja maailman talous elpyy
%% s5 = Korona kukistettu ja maailma lamassa
%% s6 = Korona kukistettu ja maailman talous elpyy
%%
%% Skenaarioiden todennakoisyydet ovat tuntemattoman koronatodennakoisyyden p0
%% funktioina
%% 
%% p1 = p0/4
%% p2 = p0/4
%% p3 = p0/4
%% p4 = p0/4
%% p5 = (1-p0)/2
%% p6 = (1-p0)/2
%%
%% Koronatodennakoisyys p0 on herkkyysparametri.
%%
%% Stefan perustaa paatoksensa odotusarvosaannolle seka kahdelle palkkiolle:
%%  
%% - Varallisuuden/lainan maara syksylla 2022
%% - Jaljella olevat opiskeluvuodet syksylla 2022
%%
%% Palkkiot oletetaan (hieman eparealistisesti) tunnetuiksi.
%%
%% Palkkioita vastaavien odotusarvosaantojen yhdistamispaino w on herkkyys-
%% parametri.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Paatosmatriisit 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Vuoden aikana kertynyt varallisuus/laina.

or = 12*252.76;      %% Opintoraha
os = 5000;           %% Opiskelijavaihtostipendi Kiinaan
km = 2000;           %% Kiinan muuttokustannukset
ms = 12*500;         %% Elamisen kulut Suomessa
mk = 12*250;         %% Elamisen kulut Kiinassa

tk = or+os-km;       %% Opiskelijan tulot Kiinassa
ts = or;             %% Opiskelijan tulot Suomessa

rl = 12*200*0.9;     %% Ruokalahetin "nettoperuspalkka"
pt = 12*2500*0.75;   %% Paatyon "nettoperuspalkka" kun vienti vetaa

R1 = [
tk-0.8*mk      tk-mk       tk-0.8*mk      tk-mk       tk-0.9*mk     tk-mk   ;   
ts-0.6*ms      ts-ms       ts-0.6*ms      ts-ms       ts-0.8*ms     ts-ms   ;    
ts+4*rl-0.6*ms ts+4*rl-ms  ts+2*rl-0.6*ms ts+2*rl-ms  ts+rl-0.8*ms  ts+rl-ms; 
0.3*pt-1.6*ms  0.7*pt-2*ms 0.3*pt-1.6*ms  0.8*pt-2*ms 0.5*pt-1.8*ms pt-2*ms ;  
];

%% Jaljella oleva valmistumisaika syksylla 2022 (tavoiteaika on yksi vuosi!!!)

R2 = [
2 2 2 2 2 2;      
2 2 1 1 1 1;      
2 2 2 2 2 2;      
3 4 3 4 3 5;      
];

R2 = -R2;        %% Maksimointimuoto!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Pienehkot taulukot p0 ja w herkkyyksille.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% p0 ja w otteistukset
p0arvot = 0:0.25:1;
warvot = 0:0.25:1;

%% Paatosongelmien vakiokentat
po1.pmat = R1;
po1.psaanto = "oa";
po2.pmat = R2;
po2.psaanto = "oa";

%% Skaalattu paatokset V1 ja V2
for i = 1:length(p0arvot)
    p = [
         p0arvot(i)/4;
         p0arvot(i)/4;
         p0arvot(i)/4; 
         p0arvot(i)/4;
         (1-p0arvot(i))/2;
         (1-p0arvot(i))/2
        ];
    po1.tn = p;
    po2.tn = p;    
    pr1 = pmatriisi(po1);
    pr2 = pmatriisi(po2);    
    V1(:,i) = pr1.parvot;
    V1(:,i) = (V1(:,i) - min(V1(:,i)))/(max(V1(:,i))-min(V1(:,i)));
    V2(:,i) = pr2.parvot;
    V2(:,i) = (V2(:,i) - min(V2(:,i)))/(max(V2(:,i))-min(V2(:,i)));
end

%% Arvojen tasapainoinen (w=0.5) yhdistaminen, kun p0=0:0.25:1
V = 0.5*V1 + 0.5*V2;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Visualisointi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Otteistustaajuus.  Mita isompi sen parempi, mutta laskeminen kestaa.
%% Jos N=500 laskeminen kestaa vajaan minuutin.
N=100;      

%% Laskut
p0arvot = linspace(0,1,N);
warvot = linspace(0,1,N);
z = zeros(N);

for i=1:length(p0arvot)
    for j=1:length(warvot)
        p0 = p0arvot(i);
        w = warvot(j);
        q0 = 1-p0;
        po1.tn = [p0/4 p0/4 p0/4 p0/4 q0/2 q0/2]';
        po2.tn = [p0/4 p0/4 p0/4 p0/4 q0/2 q0/2]';
        pr1 = pmatriisi(po1);
        pr2 = pmatriisi(po2);
        v1 = pr1.parvot;
        v1 = (v1-min(v1))/(max(v1)-min(v1));
        v2 = pr2.parvot;
        v2 = (v2-min(v2))/(max(v2)-min(v2));
        v = w*v1 + (1-w)*v2;
        [roskaa, z(j,i)] = max(v);
    end    
end
    
%% Itse kuva    
contourf(p0arvot, warvot, z,1:4)
set(gca, "fontsize", 24);
xlabel("p_0")
ylabel("w")
puna  = [1 0 0];
viher = [0 1 0];
colormap([puna; 0.6*puna; viher; 0.6*viher])
%% a1 = puna
%% a2 = tumma puna
%% a3 = viher
%% a4 = tumma viher
text(0.40,0.50,"a_1", "fontsize",24)
text(0.20,0.10,"a_2", "fontsize",24)
text(0.85,0.60,"a_3", "fontsize",24)
text(0.10,0.85,"a_4", "fontsize",24)
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Herkkyyskuvan piirto PDF-tiedostoon.  Poista kommentit, jos haluat PDF:n
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%orient("landscape");
%print("ss_herkkyyskuva.pdf");