SJÄLVSTÄNDIGA ARBETEN I MATEMATIK MATEMATISKA INSTITUTIONEN, STOCKHOLMS UNIVERSITET

SJÄLVSTÄNDIGA ARBETEN I MATEMATIK

MATEMATISKA INSTITUTIONEN, STOCKHOLMS UNIVERSITET

Modellering av omställning av odling med markovkedjor

av

Henrik Hansson

2018 - No K22

Modellering av omställning av odling med markovkedjor

Henrik Hansson

Självständigt arbete i matematik 15 högskolepoäng, grundnivå Handledare: Yishao Zhou

Sammanfattning

Detta arbete ¨ar en till¨ampning av markovkedjor f¨or modellering av f¨or¨andringar av odlingsm¨onster.

Markovkedjor ¨ar en matematisk modell som bygger p˚a att man definierar sannolik- heter f¨or f¨orflyttningar mellan olika tillst˚and. Vissa markovkedjor har egenskapen att efter tillr¨ackligt l˚ang tid uppn˚a ett station¨art tillst˚and. Den matematiska teorin f¨or station¨ara tillst˚and av markovkedjor beskrivs kortfattat i detta arbete.

Overg˚¨ angssannolikheter att odla viss gr¨oda och odlingsmetod f¨or n¨astkommande odlingss¨asong har tagits fram. Genom att r¨akna ut markovkedjor har f¨or¨andringen av odlingsmetod fr˚an traditionell till ekologisk odling kunnat studeras. En modell har tagits fram d¨ar n¨arliggande byar kan p˚averka f¨or¨andringsben¨agenheten av val av gr¨oda och odlingsmetod. Det inneb¨ar att ¨overg˚angssannolikheterna f¨or mar- kovkedjorna f¨or¨andras under simuleringen. Resultatet av simuleringarna visar att globala f¨or¨andringar av ¨overg˚angssannolikheter ger ett snabbt genomslag av val av odlingsmetod medan det tar l¨angre tid innan lokala f¨or¨andringar f˚ar e↵ekt. P˚a l¨angre sikt kan lokala stimulanser som sprider sig ge stora f¨or¨andringar. I detta fall val av gr¨oda och odlingsmetod.

Inneh˚ all

1 Inledning 3

1.1 Syfte . . . 3

1.2 Problemformulering . . . 3

2 Teori 4 2.1 Allm¨an teori om markovkedjor . . . 4

2.2 Egenv¨arden och egenvektorer . . . 7

2.3 Grundl¨aggande teori om singul¨ara matriser . . . 8

2.4 Utr¨akning av fixerade radvektorer och kolumnvektorer . . . 9

2.5 N˚agra egenskaper hos stokastiska matriser . . . 10

3 Genomf¨orande 13 3.1 Beskrivning av och m˚al f¨or modellen . . . 13

3.2 Deluppgift 1 - statiska markovkedjor . . . 16

3.3 Deluppgift 2 - dynamiska markovkedjor . . . 17

4 Resultat 19 4.1 Deluppgift 1 - statiska markovkedjor . . . 19

4.2 Deluppgift 2 - dynamiska markovkedjor . . . 22

5 Diskussion 25

6 Referenser 26

1 Inledning

1.1 Syfte

Syftet med detta arbete ¨ar att till¨ampa teorin om markovkedjor genom att skapa en matematisk modell som inneh˚aller markovkedjor. Utifr˚an modellen ska sedan ber¨akningar och simuleringar ske med det statistiska verktyget R.

1.2 Problemformulering

Detta arbete ska visa hur snabbt omst¨allningen av gr¨odor och ¨overg˚angen till eko- logisk odling sker med hj¨alp av en matematisk modell baserad p˚a markovkedjor.

Overg˚¨ angsmatriser f¨or markovkedjor ska byggas upp av tv˚a typer av parametrar.

Den ena typen av parameter ¨ar preferens f¨or viss gr¨oda eller odlingsmetod. Den andra parametern ¨ar f¨or¨andringstr¨ogheten f¨or odlingstyper. Simuleringarna ska ocks˚a visa vad som h¨ander vid lokala respektive globala insatser. E↵ekten av att n¨arliggande g˚ardar/byar p˚averkar varandra ska modelleras och simuleras.

Modellen bygger p˚a att till¨ampa markovkedjor. Diskreta markovkedjor har egen- skapen att de f¨or¨andras stegvis och ofta uppn˚ar station¨ara tillst˚and.

2 Teori

2.1 Allm¨an teori om markovkedjor

En markovkedja ¨ar en tidsdiskret stokastisk process. Markovkedjan kan anta ett visst antal tillst˚and i m¨angden S ={s¹, s2, . . . , sr}. Processen b¨orjar i ett av till- st˚anden si och f¨orflyttar sig sedan till ett annat tillst˚and sj med en sannolikhet sij. Processen kan ocks˚a stanna kvar i samma tillst˚and med sannolikheten sii. Overg˚¨ angssannolikheterna beror inte p˚a vilka tidigare tillst˚and som kedjan befun- nit sig i. ¨Overg˚angssannolikheterna kan representeras med en ¨overg˚angsmatris P.

Nedan visas ett exempel p˚a en ¨overg˚angsmatris.

0

@

s1 s2 s3

s1 0.2 0.4 0.4 s2 0.5 0.25 0.25 s3 0.3 0.7 0

1 A

L˚at p⁽ⁿ⁾_ij beteckna sannolikheten att om man startar i tillst˚and i hamnar i tillst˚and j efter n steg. D˚a ges sannolikheten att man befinner sig i ett visst tillst˚and tv˚a tillst˚and senare av:

p⁽²⁾₁₃ = p11p13+ p12p23+ p13p33

Detta uttryck f˚as ¨aven om man multiplicerar rad ett i P med kolumn tre i P.

Generellt om markovkedjan har r olika tillst˚and s˚a g¨aller:

p⁽²⁾_ij = Xr

k=1

pikpkj

Teorem 1

L˚at P vara ¨overg˚angsmatrisen f¨or en markovkedja. Elementet ij i p⁽ⁿ⁾_ij av matrisen Pⁿ ger sannolikheten att markovkedjan befinner sig i tillst˚and sj efter n steg med start i tillst˚and si.

Bevis Xr

k=1

pikpkj= p⁽²⁾_ij ¨ar element (i,j) i matrismultiplikationen PP = P² f¨or P med r

rader och kolumner.

Xr k=1

p⁽²⁾_ik pkj = p⁽³⁾_ij ¨ar element (i,j) i matrismultiplikationen

P²P = P³. Forts¨atter man s˚a h¨ar s˚a kommer man till att Xr

k=1

p^{(n 1)}_ik pkj = p⁽ⁿ⁾_ij ¨ar element (i,j) i matrismultiplikationen P^{n 1}P = Pⁿ

Teorem 2

L˚at P vara ¨overg˚angsmatrisen f¨or en markovkedja och l˚at u vara sannolikhetsvek- torn f¨or startf¨ordelningen. Sannolikheten att markovkedjan befinner sig i tillst˚and si efter n steg ges av det i:te vektorelementet i f¨oljande vektor:

u⁽ⁿ⁾ = uPⁿ Bevis

u⁽¹⁾_j = Xr k=1

ukpkj ¨ar sannolikheten att man befinner sig i tillst˚and uj efter ett steg.

Samma uttryck g¨aller f¨or vektorelement u⁽¹⁾_j i matrismuliplikationen u⁽¹⁾= uP.

u⁽²⁾_j = Xr k=1

ukp⁽²⁾_kj motsvarar vektorelement u⁽²⁾_j i matrismultiplikationen u⁽²⁾= uP².

Forts¨atter man s˚a h¨ar s˚a kommer man till att u⁽ⁿ⁾_j = Xr

k=1

ukp⁽ⁿ⁾_kj vilket motsvarar vektorelement u⁽ⁿ⁾_j i matrismultiplikationen u⁽ⁿ⁾= uPⁿ.

Exempel

Tv˚a fotbollslag i h¨ogsta serien har m¨ott varandra m˚anga g˚anger genom ˚aren. Lagen kallas f¨or Tigrarna och Vargarna. Nedan visas ¨overg˚angsmatrisen P f¨or sannolik- heten, beroende p˚a resultatet i den senaste matchen, att utfallet blir vinst f¨or Tigrarna (T), vinst f¨or Vargarna (V) eller oavgjort (O).

P = 0

@

T V O

T 0.60 0.20 0.20 V 0.20 0.30 0.50 O 0.30 0.20 0.50

1 A

Sannolikheterna f¨or hur det g˚ar mellan lagen efter tv˚a matcher (d.v.s. utan vetskap om resultatet f¨or matchen innan den aktuella matchen) ges av P², vilket visas nedan.

P²= 0

@

T V O

T 0.46 0.22 0.32 V 0.33 0.23 0.44 O 0.37 0.22 0.41

1 A

Sannolikheterna f¨or hur det g˚ar mellan lagen efter tre matcher (d.v.s. utan vetskap om resultatet f¨or de tv˚a matcherna innan den aktuella matchen) ges av P³, vilket visas nedan.

P³ = 0

@

T V O

T 0.416 0.222 0.362 V 0.376 0.223 0.401 O 0.389 0.222 0.389

1 A

Sannolikheterna f¨or hur det g˚ar mellan lagen efter n matcher ¨ar limn!1Pⁿ, vilket visas nedan.

limn!1Pⁿ = 0

@

T V O

T 0.397 0.222 0.381 V 0.397 0.222 0.381 O 0.397 0.222 0.381

1 A

Man kan notera att raderna ¨ar likadana vilket inneb¨ar att om man multiplicerar vilken startvektor u som helst med limn!1Pⁿ s˚a blir resultatet radvektorn i ma- trisen limn!1Pⁿ.

F¨or vidare diskussion i denna text beh¨over vi n˚agra beteckningar:

Matrisen P s¨ages vara icke-negativ om alla element i P ¨ar icke negativa, vilket betecknas P 0.

Matrisen P s¨ages vara positiv om alla element ¨ar positiva, vilket betecknas P > 0.

Matrisen P kallas f¨or en stokastisk matris om alla radsummor = 1.

L˚at c = [1, 1, . . . , 1]^T. Definition 1

En marovkedja s¨ags vara ergodisk om det ¨ar m¨ojligt att g˚a fr˚an varje tillst˚and till varje tillst˚and (utan att n¨odv¨andigtvis g¨ora det i ett steg).

Definiton 2

En markovkedja s¨ags vara regulj¨ar om P^m> 0 f¨or n˚agot positivt heltal m.

Med andra ord inneb¨ar det att man f¨or n˚agot n kan g˚a fr˚an vilket tillst˚and som helst till vilket tillst˚and som helst i exakt n steg. Detta visar att varje regulj¨ar mar- kovkedja ¨ar ergodisk. Men en ergodisk markovkedja ¨ar n¨odv¨andigtvis inte regulj¨ar.

En ¨overg˚angsmatris som inte har n˚agra nollor definierar en regulj¨ar markovkedja.

Men det finns ¨aven regulj¨ara markovkedjor som har nollor i ¨overg˚angsmatrisen.

Ett exempel visas nedan d˚a

P = 0

@1/2 1/2 0 1/2 1/4 1/4

0 1/2 1/2 1 A

¨ar regulj¨ar eftersom

P² = 0

@1/2 3/8 1/8 3/8 7/16 3/16 1/4 3/8 3/8

1 A

Teorem 3

F¨or en ergodisk markovkedja g¨aller att limn!1p⁽ⁿ⁾_ij = wi d¨ar w ={w¹, w2, ..., wm}

¨ar en strikt positiv sannolikhetsvektor. Bevis f¨or denna sats finns i avsnitt 2.5.

Definition 3

En radvektor w f¨or vilken wP = w kallas w f¨or en fixerad radvektor f¨or P. P˚a liknande s¨att om x ¨ar en kolumnvektor s˚a att Px = x kallas x f¨or en fixerad kolumnvektor f¨or P.

Man kan r¨akna ut fixerade rad- och kolumnvektorer f¨or en regulj¨ar markovked- ja p˚a flera olika s¨att. Ett s¨att att r¨akna fram vektorerna ¨ar att r¨akna fram h¨oger- och v¨ansteregenvektorer till ¨overg˚angsmatrisen. I f¨oljande avsnitt beskrivs bl.a.

teorin f¨or detta.

2.2 Egenv¨arden och egenvektorer

Egenvektorer och egenv¨arden till en matris A definieras genom icke triviala l¨osningar (x 6= 0) till ekvationen Ax = x. Detta kan skrivas som [A I]x = 0 d¨ar I ¨ar identitetsmatrisen. L¨osningen till det(A I) = 0 ger det karakteristiska polynomet f¨or vilka egenv¨ardena kan r¨aknas ut. Egenvektorer (korrespondera- de till egenv¨ardet) r¨aknas sedan ut genom att s¨atta in egenv¨ardet i relationen

A I = 0. N¨ar x st˚ar till h¨oger om A kallas egenvektorerna och egenv¨ardena f¨or h¨ogeregenvektorer och h¨ogeregenv¨arden h. N¨ar relationen ist¨allet ser ut s˚a h¨ar vA = v s˚a kallas egenvektorerna och egenv¨ardena f¨or v¨ansteregenvektorer och v¨ansteregenv¨arden v.

L˚at A vara en n x n matris med element aij. L˚at vidare Ri(A) =|P

i6=jaij|. D˚a ¨ar Gershgorins diskskiva den diskskiva D(aii, Ri(A)) centrerad i aii p˚a det komplexa planet med radie Ri(A).

Teorem 4

Varje egenv¨arde av en matris ligger i ˚atminstone en Gershgorins diskskiva. Det vill s¨aga om en n x n matris A har ett egenv¨arde , d˚a g¨aller att:

2 [n i=1

D(aii, Ri(A))

Beviset f¨or denna sats h¨anvisas till ref.[2].

Exempel 1

P ¨ar en stokastisk matris och pij 0 D˚a g¨aller Ri(P ) =P

i6=jpij = 1 pii. S˚a alla egenv¨arden ligger i Sn

i=1D(pii, 1 pii).

Mer precist g¨aller enligt omv¨anda triangelolikheten att

|z| pii  |z pii|  1 pii, |z|  1 ) 1 ¨ar det st¨orsta egenv¨ardet till P Exempel 2

P = P I d¨ar > 0 s˚a att alla element i P I ¨ar positiva.

Nu ¨ar Ri(P ) = Ri(P ). Gershgorins diskskivor ¨ar

|z (pii )|  Ri(P ) = 1 pii

) |z| (pii ) 1 pii ) |z|  1

S˚a alla egenv¨arden av P ¨ar begr¨ansade av 1 som ¨ar ett egenv¨arde till P .

2.3 Grundl¨aggande teori om singul¨ara matriser

Definition 4

En matris A ¨ar singul¨ar om det existerar en vektor v 6= 0 s˚a att Av = 0.

Teorem 5

¨ar ett egenv¨arde till P om och endast om P I ¨ar singul¨ar f¨or n˚agon v6= 0.

Bevis

Definitionen av egenv¨arde ¨ar att det existerar en skal¨ar s˚adan att P x = x ) (P I)x = 0. F¨or att det ska finnas en icke trivial l¨osning (x 6= 0) till denna ekvation kr¨avs enligt definitionen ovan att P I ¨ar singul¨ar.

Teorem 6

A ¨ar singul¨ar om och endast om kolumnerna i A ¨ar linj¨art beroende och om och endast om determinanten det(A) = 0.

Bevis

Bevis f¨or denna sats h¨anvisas till ref.[2].

2.4 Utr¨akning av fixerade radvektorer och kolumnvektorer

F¨or alla stokastiska ¨overg˚angsmatriser g¨aller att P

jPij = 1) P c = c. Detta vi- sar att c ¨ar en h¨ogeregenvektor till P . c multiplicerad med en godtycklig konstant k ¨ar ocks˚a en l¨osning p˚a matrisekvationen. Men eftersom vi ¨ar intresserade av en stokastisk vektor (sannolikhetsf¨ordelning) m˚aste absolutbeloppet av den s¨okta vektorn |x| = 1 ) x = c

P

jcj

. S˚a x ¨ar allts˚a den station¨ara h¨ogeregenvektorn till

¨overg˚angsmatrisen P . Den station¨ara h¨ogeregenvektorn x ¨ar den fixerade kolumn- vektorn.

x ¨ar en h¨ogeregenvektor men det vi speciellt s¨oker ¨ar v¨ansteregenvektorn v som uppfyller vP = v. Relationen f¨or v¨ansteregenv¨arden lyder vP = vv men ef- tersom det ¨ar ett station¨art tillst˚and s˚a ¨ar v = 1. Det kan d˚a skrivas som (vP )^T = v^T ) P^Tv^T = v^T ) (P^T I)v^T = 0.

0 = det(P^T I) = det(P^T I^T) = det(P I)^T = det(P I). Det sista steget eftersom det(A) = det(A^T).

Vi ser h¨ar att det karakteristiska polynomet som ges av det(P I) = 0 ¨ar likadant f¨or h¨oger- och v¨ansteregenv¨arden. D¨arf¨or g¨aller alltid att v = h. Ef- tersom det vi s¨oker ¨ar en stokastisk vektor m˚aste vi normera den s¨okta vektorn,

) w = v

P

jvj

. w ¨ar allts˚a den s¨okta station¨ara v¨ansteregenvektorn till P . Den station¨ara v¨ansteregenvektorn w ¨ar den fixerade radvektorn.

Det ¨ar m¨ojligt att P inneh˚aller nollor och ¨and˚a konvergerar till ett station¨art tillst˚and med en strikt positiv sannolikhetsvektor w (kravet ¨ar att markovkedjan

¨ar ergodisk).

2.5 N˚ agra egenskaper hos stokastiska matriser

Prop 1

Om P och P⁰ ¨ar stokastiska matriser s˚a ¨ar P P⁰ en stokastisk matris.

Bevis

L˚at P = (Pij), P⁰ = (P_ij⁰) och c = [1, 1, . . . , 1]^T. Om P och P⁰ ¨ar stokastiska d˚a ¨ar P c = c och P⁰c = c. (P P⁰)c = P (P⁰c) = P c = c. Allts˚a ¨ar P P⁰ stokastisk.

Prop 2

Om P ¨ar en stokastisk matris, s˚a har P ett egenv¨arde = 1.

Bevis

Det f¨oljer av att P ¨ar stokastisk, d.v.s. P c = c.

Fr˚an beviset ser vi att c ¨ar en egenvektor f¨or P associerad till egenv¨ardet 1.

F¨oljdsats

Om P ¨ar en stokastisk matris s˚a ¨ar Pⁿ en stokastisk matris.

Prop 3

En stokastisk matris kan inte ha egenv¨arde d¨ar | | > 1 Bevis

Se exempel 1 till teorem 4.

Vi antar nu att P har alla element positiva, d.v.s. pij > 0 f¨or alla i, j = 1, ..., n.

D˚a har vi f¨oljande teorem.

Teorem 7

Antag att P ¨ar en n x n stokastisk matris med positiva element.

(i) = 1 ¨ar ett egenv¨arde och c ¨ar en egenvektor.

(ii) = 1 ¨ar det enda egenv¨ardet p˚a enhetscirkeln.

(iii) = 1 har geometrisk multiplicitet = 1.

(iv) = 1 har algebraisk multiplicitet = 1.

Bevis

(i) ¨ar redan visat i prop 2.

(ii) Enligt prop 3, om ¨ar ett egenv¨arde till P g¨aller det att | |  1. Vi ska nu visa att olikheten ¨ar strikt. Betrakta matrisen P = P I d¨ar > 0 ¨ar tillr¨ackligt litet s˚a att P fortfarande har alla element positiva.

s˚a P c = (P I)c = P c c = c c = (1 )c ) (1 ) ¨ar ett egenv¨arde till P

L˚at y vara en egenvektor associerad till , d.v.s. P y = y ( 6= 1) P y = (P I)y = P y y = y y = ( )y

) ( ) ocks˚a ¨ar ett egenv¨arde till P

Enligt exempel 2 efter teorem 4 ¨ar 1 det st¨orsta egenv¨ardet till P . S˚a | |  1 . Nu ska vi visa att det enda egenv¨ardet till P med beloppet| | = 1 ¨ar = 1.

Det ¨ar klart att = 1 inte ¨ar ett egenv¨arde ty| 1 | = 1 + > 1 . Nu antar vi att | | = 1 men ska inte vara 1, d.v.s. = a + bi, b6= 0 och a²+ b² = 1.

) | |² = |a + bi |² = (a )²+ b² = a² 2a + ²+ b² = 1 2a + ² >

1 2 + ² = (1 )²

) | | > 1 , en mots¨agelse. Det finns inte andra egenv¨arden p˚a enhetscirkeln.

(iii) Antag geometrisk multiplicitet k > 1 av ett egenv¨arde = 1, d.v.s. dimen- sionen till nollrummet av matrisen I P ¨ar st¨orre ¨an 1. Med andra ord finns en annan egenvektor y 6= c som ¨ar linj¨art oberoende av c, associerad till egenv¨ardet 1 : P y = y. y m˚aste vara en reell vektor eftersom P ¨ar en reell matris. Vi har c > 0 (komponetvis). D˚a ¨ar det m¨ojligt att hitta en w = ↵c + y f¨or n˚agot ↵2 R s˚adan att n˚agon komponent i w = 0. P w = ↵P c + P y = ↵c + y = 1w = w.

Men P > 0 s˚a P w > 0 vilket inneb¨ar att w > 0. Det ¨ar en mots¨agelse eftersom vi antog att w inte ¨ar strikt positiv. Allts˚a kan = 1 inte ha tv˚a linj¨art oberoende egenvektorer, d.v.s. den geometriska multipliciteten m˚aste vara 1.

(iv) Fr˚an (iii) vet vi att = 1 bara har en egenvektor c > 0. Antag att c har algebraisk mulitiplicitet st¨orre ¨an 1. F¨or enkelhets skull antar vi att den algebra- iska multipliciteten ¨ar 2, d˚a ser Jordanblocket ut s˚a h¨ar

✓1 1 0 1

◆ .

D˚a finns en vektor z s˚a att (P I)z = y och (P I)y = 0. y m˚aste (eftersom P ¨ar stokastisk) vara en multipel av c s˚a vi antar att y = c utan att g¨ora n˚agon inskr¨ankning. F¨or egenv¨ardet 1 kan vi hitta en positiv vektor f s˚a att P^Tf = f (d.v.s. v¨ansteregnenvektor). D˚a f^T(P I)z = f^TP z f^Tz = f^Tz f^Tz = 0.

Men f^T(P I)z = f^Tc > 0 (eftersom f > 0, c > 0). Det inneb¨ar att algebraisk multiplicitet inte kan vara st¨orre ¨an 1.

Bevis av teorem 3

Om matrisen P ¨ar diagonaliserbar s˚a g¨aller Pⁿu = a1Pⁿu1+a2Pⁿu2+...+amPⁿum, d¨ar u1, u2, ..., um ¨ar egenvektorer och a1, a2, ..., am ¨ar konstanter.

Puk = kuk, u1= c, 1= 1) Pⁿu = a1c + a2 n

2u2+ ... + am n mum. Eftersom alla egenv¨arden 2...m< 0 ) limⁿ!1Pⁿu = a1c eller Pc = c.

Nu antar vi att matrisen P inte ¨ar diagonaliserbar. Antag P har egenv¨arden 1, 2, ... k d¨ar j < 1 f¨or j = 2, ..., k. S˚a P har Jordans kanonisk form J =

diag (J1, ...Jk) d¨ar J1= 1 och Jj f¨or j = 2, .., k ser ut s˚a h¨ar:

Jj = 0 BB B@

j 1

. .. ...

. .. 1

j

1 CC

CA, j = 2, ..., k

Observera att storlek p˚a dessa block beror p˚a den geometriska multipliciteten f¨or

j. Att bevisa Pⁿu konvergerar ¨ar samma som att bevisa Jⁿv konvergerar. Men Jⁿ = diag (J₁ⁿ, ...J_kⁿ)! diag (1, 0, ...0, ) eftersom

J_jⁿ = 0 BB B@

nj n 1

n 1

j · · · _{l 1}^{n n l+1}j

. .. . ..

. .. ⁿ₁ ^{n 1}_j

nj

1 CC

CA, j = 2, ..., k

och | j| < 1. H¨ar har vi antagit att blockstorlek f¨or Jj ¨ar l. D¨armed ¨ar det klart med beviset f¨or teorem 3.

3 Genomf¨orande

3.1 Beskrivning av och m˚ al f¨or modellen

Modellen bygger p˚a att anv¨anda markovkedjor. Diskreta markovkedjor har egen- skapen att f¨or¨andringar sker stegvis och ofta uppn˚as station¨ara tillst˚and. Denna modell bygger p˚a att definiera en ¨overg˚angsmatris P f¨or varje bonde, by, stad eller region. Varje ¨overg˚angsmatris ˚aterspeglar den lokala situationen f¨or den enskilda bonden, byn, staden eller regionen.

Vi kan t¨anka oss att det g¨aller Kina och att vi har tre odlingss¨asonger per ˚ar. Det odlas tre olika s¨adesslag, ris, vete och kik¨artor. Det ¨ar milj¨o- och klimatm¨assigt f¨ordelaktigt att odla de olika gr¨odorna i f¨oljande ordning: kik¨artor, vete och ris.

Varje gr¨oda kan odlas traditionellt och ekologiskt. Varje odlingss¨asong represente- rar ett steg i markovkedjan.

En 6x6 stor ¨overg˚angsmatris P definieras. De sex olika tillst˚anden (s1 s6) i P av- ser vilken odlingstyp (OT) som odlas. Odlingstypen anges med tv˚a bokst¨aver. Den f¨orsta bokstaven anger gr¨oda (K=kik¨artor, V=vete, R=ris) och den andra bok- staven anger odlingsmetod (T=traditionell, E=ekologisk). Tex. betyder KT tradi- tionell odling av kik¨artor. Tillst˚anden definieras s˚a att s1=KT, s2=KE, s3=VT, s4=VE, s5=RT och s6=RE. Matriselementen i P ska tolkas som sannolikheten pij att en odlingstyp ¨overg˚ar till en annan odlingstyp eller forts¨atter med samma odlingstyp efter n¨asta sk¨ord (odlingss¨asong). Vilken geografisk spridning som en odlingstyp avser i ¨overg˚angsmatrisen kan variera beroende p˚a vad man vill simu- lera. Exempel ¨ar odling p˚a samma jordbruksmark, odling av samma bonde eller odling i samma by eller stad.

Gr¨odornas preferenser definieras genom en viktning s˚a att summan av gr¨odornas vikter = 100%. Ex. vikt kik¨artor = 20%, vikt vete = 30% och vikt ris = 50%.

Varje gr¨oda viktas sedan p˚a odlingsmetod (ekologisk eller traditionell odling). Ex.

vikt trad = 80%, vikt eko = 20%. De tv˚a vikterna multipliceras med varandra f¨or att f˚a f¨ordelningen f¨or de 6 olika tillst˚anden. Ex. vikt ris trad = 0,5*0,8=0,4.

Det finns en tr¨oghet i odlingen som beror p˚a traditioner som g¨or att en jordbruka- re g¨arna forts¨atter att odla samma gr¨oda med samma odlingsmetod (traditionellt eller ekologiskt). Tr¨ogheten T definieras som en faktor i intervallet [0 1]. 1 = max- imal tr¨oghet och 0 = ingen tr¨oghet. Den del som ska f¨ordelas mellan gr¨odorna enligt viktningen r¨aknas ut genom att multiplicera med 1-T. Icke diagonalelement

i P f˚ar v¨arden enligt denna f¨ordelning. F¨or diagonalelementen adderas T.

Nedan visas ett exempel p˚a en ¨overg˚angsmatris P d¨ar tr¨ogheten T = 0,9.

0 BB BB BB

@

KT KE V T V E RT RE

KT 0, 916 0, 004 0, 024 0, 006 0, 04 0, 01 KE 0, 016 0, 904 0, 024 0, 006 0, 04 0, 01 V T 0, 016 0, 004 0, 924 0, 006 0, 04 0, 01 V E 0, 016 0, 004 0, 024 0, 906 0, 04 0, 01 RT 0, 016 0, 004 0, 024 0, 006 0, 94 0, 01 RE 0, 016 0, 004 0, 024 0, 006 0, 04 0, 91

1 CC CC CC A

Rad fem i P r¨aknas fram enlig f¨oljande:

p51 = vikt kik¨artor * vikt trad * (1-T) = 0,2*0,8*0,1 = 0,016 p52 = vikt kik¨artor * vikt eko * (1-T) = 0,2*0,2*0,1 = 0,004 p53 = vikt vete * vikt trad * (1-T) = 0,3*0,8*0,1 = 0,024 p54 = vikt vete * vikt eko * (1-T) = 0,3*0,2*0,1 = 0,006

p55 = vikt ris * vikt trad * (1-T) + T = 0,5*0,8*0,1 + 0,9 = 0,94 p56 = vikt ris * vikt eko * (1-T) = 0,5*0,2*0,1 = 0,01

N¨ar en ¨overg˚angsmatris skapas p˚a detta s¨att kallar vi den grundmatris Pg. Fr˚an Pg kan man g¨ora f¨or¨andringar f¨or specifika ¨overg˚angssannolikheter f¨or att se hur dessa f¨or¨andringar p˚averkar simuleringsresultatet.

Simuleringarna ska visa hur snabbt odlingstypen f¨or¨andras vid olika f¨or¨andringar av de tv˚a olika parametertyperna preferenser och tr¨oghet. Vilka faktorer p˚averkar preferensen av odlingstyp? N˚agra exempel ¨ar pris p˚a t.ex. uts¨ade, konstg¨odsel, bek¨ampningsmedel, efterfr˚agan, jordm˚an, teknisk utveckling, milj¨okonsekvenser (t.ex. ¨overg¨odda vattendrag). Det finns m˚anga s¨att att f¨or¨andra preferenserna f¨or en odlingstyp, t.ex. f¨orbjuda vissa bek¨ampningsmedel eller h¨oja skatter.

Vilka faktorer p˚averkar tr¨ogheten f¨or en odlingstyp? N˚agra exempel ¨ar traditio- ner, utbildningsniv˚a, kunskap och omst¨allningskostnader. Hur kan man p˚averka tr¨ogheten? Ett exempel ¨ar att ge subventioner till odling av viss gr¨oda eller sub- ventioner f¨or omst¨allning fr˚an konventionellt jordbruk till ekologiskt.

Den f¨orsta deluppgiften ¨ar att simulera f¨or¨andringarna utan att n¨arliggande by- ar (markovkedjor) p˚averkar varandra. Det station¨ara tillst˚andet ska r¨aknas fram och fixerade radvektorer ska r¨aknas ut. Simuleringarna ska visa vilka skillna-

derna blir n¨ar man f¨or¨andrar preferenser (t.ex. ¨okad efterfr˚agan) j¨amf¨ort med n¨ar man f¨or¨andrar tr¨ogheten (t.ex. h¨ojd subvention). Olika starttillst˚and (initiala f¨ordelningen) ska simuleras. Tiden det tar (med viss noggrannhet) tills det sta- tion¨ara tillst˚andet uppn˚atts ska simuleras. F¨or¨andringstakten ska r¨aknas ut.

I den andra deluppgiften ska intilliggande markovkedjor p˚averka varandra. En in- tressant fr˚aga att simulera ¨ar skillnaderna mellan global respektive lokal p˚averkan.

T.ex. hur p˚averkar mindre globala subventioner j¨amf¨ort med st¨orre lokala sub- ventioner? Grannb¨onder och n¨arliggande byar p˚averkar varandra. Traditioner kan brytas och kunskap kan ¨overf¨oras. F¨or att kunna simulera detta inf¨ors smittoef- fekter i den andra delen av simuleringen.

I denna modell har varje markovkedja (g˚ard , by eller stad) fyra grannar (ej di- agonala grannar) i ett regelbundet rutn¨at. Smittoe↵ekter ¨overf¨ors till samtliga grannar. Smittoe↵ekten definieras p˚a s˚a s¨att att tr¨ogheten minskar en viss procent (smittofaktorn) p˚a skillnaden mellan tr¨ogheten f¨or tv˚a n¨arliggande g˚ardar. Smit- toe↵ekter r¨aknas f¨or varje granne och adderas upp. I dessa simuleringar kommer tr¨ogheter aldrig ¨oka p˚a grund av smittoe↵ekter. Ex. en g˚ard (med tr¨oghet 0,96) har fyra grannar med l¨agre tr¨oghet (0.92, 0.90, 0.94 och 0.88). Smittofaktor = 5%

) nya tr¨ogheten f¨or g˚arden = 0,96 - 0,05*(0,04 + 0,06 + 0,02 + 0,08) = 0,95.

Smittober¨akningar g¨ors inf¨or varje ny sk¨ord (eller efter ett visst antal sk¨ordar).

Simuleringarna ska visa vilka skillnader som uppst˚ar n¨ar man f¨or¨andrar tr¨ogheten globalt j¨amf¨ort med n¨ar man f¨or¨andrar tr¨ogheten lokalt med smittoe↵ekter.

Simuleringar som ska j¨amf¨oras med varandra ska vara ekonomiskt likv¨ardiga, t.ex.

den totala subventionen. P˚a s˚a s¨att kan den samh¨allsekonomiska nyttan studeras f¨or olika ˚atg¨arder.

Dessa simuleringar har ingen verklighetsanknytning eftersom modellen ¨ar allde- les f¨or enkel och all data ¨ar helt och h˚allet fiktiv. Det simuleringarna kan ge ¨ar en uppfattning om hur olika typer av f¨or¨andringar i modellen p˚averkar den l˚angsiktiga f¨or¨andringen i odlingsm¨onster.

Modelleringen, ber¨akningarna och simuleringarna kommer att ske i R.

3.2 Deluppgift 1 - statiska markovkedjor

Resultatet av ber¨akningarna redovisas f¨or fem olika (representativa) ¨overg˚angsmatriser.

Valet av matriser har gjorts med syftet att visa f¨or¨andringarna n¨ar man ¨andrar

¨overg˚angssannolikheter mellan odlingsmetoderna (traditionellt och ekologiskt) och n¨ar man f¨or¨andrar tr¨ogheten. F¨or¨andringar i ¨overg˚angssannolikheter mellan od- lingstyper med samma odlingsmetod visade inga avvikande tendenser s˚a d¨arf¨or redovisas inte dessa resultat. Att det f¨orh˚aller sig s˚a ¨ar fullt logiskt eftersom inga s¨arskilda modellegenskaper tagits fram f¨or ¨overg˚angssannlikheter mellan gr¨odor.

P1=

0 BB BB BB

@

KT KE V T V E RT RE

KT 0, 916 0, 004 0, 024 0, 006 0, 04 0, 01 KE 0, 016 0, 904 0, 024 0, 006 0, 04 0, 01 V T 0, 016 0, 004 0, 924 0, 006 0, 04 0, 01 V E 0, 016 0, 004 0, 024 0, 906 0, 04 0, 01 RT 0, 016 0, 004 0, 024 0, 006 0, 94 0, 01 RE 0, 016 0, 004 0, 024 0, 006 0, 04 0, 91

1 CC CC CC A

Den f¨orsta P1 ¨ar identisk med grundmatrisen Pg.

P2=

0 BB BB BB

@

KT KE V T V E RT RE

KT 0, 748 0, 012 0, 072 0, 018 0, 12 0, 03 KE 0, 048 0, 712 0, 072 0, 018 0, 12 0, 03 V T 0, 048 0, 012 0, 772 0, 018 0, 12 0, 03 V E 0, 048 0, 012 0, 072 0, 718 0, 12 0, 03 RT 0, 048 0, 012 0, 072 0, 018 0, 82 0, 03 RE 0, 048 0, 012 0, 072 0, 018 0, 12 0, 73

1 CC CC CC A

P2 ¨ar lika med Pg f¨orutom att tr¨ogheten = 0,7.

P3=

0 BB BB BB

@

KT KE V T V E RT RE

KT 0, 916 0, 004 0, 024 0, 006 0, 04 0, 01 KE 0, 012 0, 908 0, 018 0, 012 0, 03 0, 02 V T 0, 016 0, 004 0, 924 0, 006 0, 04 0, 01 V E 0, 012 0, 008 0, 018 0, 912 0, 03 0, 02 RT 0, 016 0, 004 0, 024 0, 006 0, 94 0, 01 RE 0, 012 0, 008 0, 018 0, 012 0, 03 0, 92

1 CC CC CC A

I P3har ¨overg˚angssannolikheterna mellan ekologiska odlingstyper (p22, p24, p26, p42, p44, p46, p62, p64, och p66) dubblerats i f¨orh˚allande till Pg. ¨Overg˚angssannolikheterna fr˚an ekologiska odlingstyper till traditionella har minskat proportionerligt. Tr¨ogheten

= 0,9.

P4=

0 BB BB BB

@

KT KE V T V E RT RE

KT 0, 916 0, 004 0, 024 0, 006 0, 04 0, 01 KE 0, 042 0, 916 0, 006 0, 024 0, 01 0, 04 V T 0, 016 0, 004 0, 924 0, 006 0, 04 0, 01 V E 0, 042 0, 016 0, 006 0, 924 0, 01 0, 04 RT 0, 016 0, 004 0, 024 0, 006 0, 94 0, 01 RE 0, 004 0, 016 0, 006 0, 024 0, 01 0, 94

1 CC CC CC A

I P4har ¨overg˚angssannolikheterna mellan ekologiska odlingstyper (p22, p24, p26, p42, p44, p46, p62, p64, och p66) fyrdubblats i f¨orh˚allande till Pg. ¨Overg˚angssannolikheterna fr˚an ekologiska odlingstyper till traditionella har minskat proportionerligt. Tr¨ogheten

= 0,9.

P5=

0 BB BB BB

@

KT KE V T V E RT RE

KT 0, 832 0, 008 0, 048 0, 012 0, 08 0, 02 KE 0, 024 0, 816 0, 036 0, 024 0, 06 0, 04 V T 0, 032 0, 008 0, 848 0, 012 0, 08 0, 02 V E 0, 024 0, 016 0, 036 0, 824 0, 06 0, 04 RT 0, 032 0, 008 0, 048 0, 012 0, 88 0, 02 RE 0, 024 0, 016 0, 036 0, 024 0, 06 0, 84

1 CC CC CC A

I P5har ¨overg˚angssannolikheterna mellan ekologiska odlingstyper (p22, p24, p26, p42, p44, p46, p62, p64, och p66) dubblerats i f¨orh˚allande till Pg. ¨Overg˚angssannolikheterna fr˚an ekologiska odlingstyper till traditionella har minskat proportionerligt. Tr¨ogheten

= 0,8.

F¨or samtliga fem ¨overg˚angsmatriser har P^k_i och uP^k_i r¨aknats fram f¨or k ={1, 2, ..., n}.

n har valts s˚a pass stort att markovkedjan n¨armat sig det station¨ara tillst˚andet.

Olika startvektorer u har testats men den startvektor som ¨ar mest intressant att studera och som ger den st¨orsta f¨or¨andringstakten ¨ar n¨ar man startar med 100%

traditionell odling. utrad = (0.2, 0, 0.3, 0, 0.5, 0) ¨ar den startvektor som anv¨ants mest i uppgift 1.

Station¨ara tillst˚and och fixerade radvektorer wi har r¨aknats fram.

Metoden som anv¨ands i R f¨or att r¨akna ut de fixerade radvektorerna ¨ar att r¨akna fram egenvektorer till ¨overg˚angsmatriserna enligt avsnitt 2.4.

3.3 Deluppgift 2 - dynamiska markovkedjor

Jag har utg˚att fr˚an grundmatrisen Pg i denna uppgift. ¨Overg˚angsmatriserna har skapats genom att f¨or¨andra tr¨ogheten hos Pg. En sv˚arighet med denna modellering

¨ar hur tr¨ogheten ska f¨oras ¨over mellan ¨overg˚angsmatriserna och hur tr¨ogheten ska f¨or¨andras vid p˚averkan fr˚an n¨arliggande byar. Mina simuleringar visade att man l¨att hamnade mellan tv˚a ytterligheter. Den ena ytterligheten var att p˚averkan fr˚an n¨arliggande byar snabbt planade ut och d¨armed blev smittoe↵ekterna mel- lan byarna ringa. Den andra ytterligheten var att den f¨or¨andrade tr¨ogheten blev s˚a stor att konvergensen hos vissa ¨overg˚angsmatriser skenade exponentiellt och d¨armed kollapsade modellen. Det jag fick g¨ora var att anpassa modellparametrar f¨or att f˚a fram resultat som ligger i linje med syftet av modellen. En modellpa- rameter som p˚averkat resultatet av modelleringen ¨ar hur diagonalelementen ska r¨aknas ut efter att simuleringen startat. N¨ar jag f¨ors¨okt anv¨anda samma metod som vid initieringen av ¨overg˚angsmatriserna uppstod urartade konvergenser, vil- ket beror p˚a att diagonalelementen blev alltf¨or sm˚a efter ett tag in i simuleringen.

Jag fick ist¨allet anv¨anda en annan metod f¨or att r¨akna fram de nya diagonalele- menten, se pseudokoden nedan. Med denna utr¨akning gav simuleringarna stabila konvergenser. Den andra sv˚arigheten var att best¨amma ett skal¨art v¨arde som de- finierar tr¨ogheten efter odlingss¨asonger > 1. Det finns flera s¨att att g¨ora detta p˚a, t.ex. medelv¨ardet av diagonalelementen i P. Att anv¨anda sig av det n¨ast st¨orsta egenv¨ardet i ¨overg˚angsmatrisen ¨ar en annan id´e. (Det n¨ast st¨orsta egenv¨ardet i

¨overg˚angsmatrisen best¨ammer konvergenshastigheten.) Det jag gjorde var att pro- va olika utr¨akningar tills jag fann utr¨akningar som st¨ammer med det man f¨orv¨antar sig av modellen. Principen f¨or utr¨akningarna visas i psudokoden nedan.

Pseudokod f¨or modellalgoritmen:

1) F¨or varje by j¨amf¨ors T (tr¨ogheten) med varje r¨atvinkligt intilliggande by.

Om T ¨ar l¨agre f¨or en intilliggande by s˚a adderas skillnaden, multiplicerat med k1, till T .

2) F¨or varje by f¨or varje rad i i P: Pi= Pi (1+T ) 3) F¨or varje by f¨or varje rad i i P: Pii= Pii P

jPij+ 1 4) F¨or varje by: T = T (T k2)

5) F¨or varje by: u^(s)= u^{(s 1)}P

6) F¨or varje by: andel eko = (u^(s)2 + u^(s)₄ + u^(s)₆ )/P

ju^(s)_j

7) G˚a till 1) tills u^(s)= u⁽ⁿ⁾ (n ¨ar maximalt antal simulerade odlingss¨asonger) k1 och k2 ¨ar parametrar som p˚averkar modellen. k1 p˚averkar hur stor smitto- e↵ekten blir mellan byar. k2 p˚averkar hur tr¨ogheten f¨or¨andras p˚a grund av smit- toe↵ekter.

4 Resultat

4.1 Deluppgift 1 - statiska markovkedjor

De fixerade radvektorerna w1 w5 f¨or respektive ¨overg˚angsmatris P1 P5 visas nedan.

w1= (0.16, 0.04, 0.24, 0.06, 0.4, 0.1) w2= (0.16, 0.04, 0.24, 0.06, 0.4, 0.1) w3 = (0.15, 0.05, 0.225, 0.075, 0.375, 0.125)

w4= (0.1, 0.1, 0.15, 0.15, 0.25, 0.25) w5 = (0.15, 0.05, 0.225, 0.075, 0.375, 0.125)

Nedan visas en plot f¨or de fem markovkedjorna P1 P5 och med startvektorn utrad. Matematiskt ¨ar det som visas i plotten utradP^k_i(0, 1, 0, 1, 0, 1)^T, vilket bety- der andelen ekologisk odling (KE, VE och RE) n¨ar markovkedjan befinner sig i steg k (= odlingss¨asong).

0 10 20 30 40 50 60

0.00.10.20.30.40.5

Andel ekologisk odling efter odlingssäsonger

odlingssäsong

andel ekologisk odling

P1 EKO=20% T=0,9 P2 EKO=20% T=0,7 P3 EKO=25% T=0,9 P4 EKO=50% T=0,9 P5 EKO=25% T=0,8

0 5 10 15 20 25 30

0.000.050.100.150.200.250.30

Andel ekologisk odling efter odlingssäsonger

odlingssäsong

andel ekologisk odling

P1 EKO=20% T=0,9 P2 EKO=20% T=0,7 P3 EKO=25% T=0,9 P4 EKO=50% T=0,9 P5 EKO=25% T=0,8

Figurerna visar (den andra figuren inzoomad) att f¨or dessa val av ¨overg˚angssannolikheter s˚a har tr¨ogheten st¨orst betydelse f¨or hur snabbt omst¨allning till ekologisk odling sker.

4.2 Deluppgift 2 - dynamiska markovkedjor

Jag gjorde m˚anga simuleringar med olika val av ing˚angsv¨arden och modellparame- terar. Den simulering som jag visar resultatet f¨or i nedanst˚aende figur har f¨oljande initiala v¨arden. 3x5 = totalt 15 byar. Alla byar startar med tr¨ogheten = 0,95 f¨orutom en by i ena h¨ornet som startar med tr¨oghet = 0,8. I figuren visas ocks˚a en kurva f¨or en konstant ¨overg˚angsmatris (som allts˚a inte p˚averkas av smittoe↵ekter) med tr¨oghet = 0,92.

0 5 10 15 20 25 30 35

0.000.050.100.150.200.250.30

Andel ekologisk odling efter odlingssäsonger

odlingssäsong

andel ekologisk odling

Fix T=0,92 T=0,8

Start T=0,95 Start T=0,95 Start T=0,95

Det man kan s¨aga om resultatet av denna simulering ¨ar att konvergenskurvorna ser ut som man kan f¨orv¨anta sig utifr˚an den modell simuleringen baserar sig p˚a.

Eftersom modellen inte ¨ar tillr¨ackligt detaljerad och att den inte baseras p˚a verkli- ga f¨orh˚allanden ¨ar det inte s˚a v¨aldigt mycket man kan utl¨asa fr˚an denna simulering.

Det som dock visar sig i denna simulering och som jag tror ocks˚a g¨aller i verklig- heten ¨ar att spridningse↵ekter tar tid innan de ger e↵ekt. Vill man f˚a en snabb f¨or¨andring s˚a ¨ar det globala subventioner som ger st¨orst e↵ekt. F¨or detta fall ¨ar det e↵ektivast att p˚averka tr¨ogheten lika mycket f¨or alla byar f¨or att f˚a en snabb omst¨allning till ekologisk odling (se kurvan f¨or konstant tr¨oghet). F¨or lokala sub- ventioner ser man att n¨ar v¨al smittoe↵ekterna har spridit sig och f˚att tid att verka s˚a kommer by efter by ikapp och f¨orbi kurvan f¨or den globala konstanta

¨overg˚angsmatrisen.

5 Diskussion

Det finns en gedigen matematisk teori kring markovkedjor och egenv¨arden. Det ¨ar enkelt att anv¨anda markovkedjor i en modell och det g˚ar fort att f˚a fram resultat genom att anv¨anda markovkedjor i simuleringar. Dessa saker talar f¨or att anv¨anda markovkedjor vid modellering. Men den grundl¨aggande fr˚agan man beh¨over st¨alla sig ¨ar om markovkedjor ¨ar l¨ampliga att anv¨anda f¨or det man vill modellera.

Markovkedjor ¨ar per definition statiska vilket g¨or att om man vill prediktera n˚agot l˚angt i framtiden s˚a kan man f˚a stora fel om ¨overg˚angssannolikheterna f¨or¨andras under tiden. (Att det ¨ar sv˚arare att prediktera l˚angt in i framtiden g¨aller f¨orvisso alla modeller.) I uppgift 2 kan ¨overg˚angssannolikheterna ¨andra sig under simuleringen vilket inneb¨ar att markovkedjorna inte uppfyller definitionen av markovkedjor. Fr˚agan ¨ar om man kan tala om markovkedjor ¨overhuvudtaget n¨ar man f¨or¨andrar ¨overg˚angssannolikheterna, men det ¨ar kanske mer av en akade- misk fr˚agest¨allning. Emellertid, att kunna f¨or¨andra ¨overg˚angssannolikheter under simuleringar ¨ar n¨odv¨andigt f¨or m˚anga f¨oreteelser man vill modellera. Man kan t¨anka sig att man b¨orjar med att anv¨anda statiska markovkedjor f¨or att skapa en enkel modell och att man sedan ut¨okar och f¨orfinar modellen. Steget fr˚an uppgift 1 till uppgift 2 ¨ar ett enkelt exempel p˚a hur det kan g˚a till.

F¨or att modellera och simulera den typ av till¨ampning som detta arbete hand- lar om (omst¨allning av odlingsmetod) p˚a riktigt s˚a tror jag man beh¨over ta fram mycket statistik och historiskt data. Man kan d˚a testa modellen mot historiskt data och se om modellen verkar ¨overrensst¨amma med verkligheten.

6 Referenser

[1] William J. Stewart , Introduction to the numeric solution of Markov chains, December 1994

[2] A. Holst och V. Ufnarovski, Matrix Theory, Studentlitteratur, 2014

[3] St´ephane Bellon, Servane Penvern, Organic Farming, Prototype for Sustainable Argicultures, 2014