{VERSION 6 0 "IBM INTEL NT" "6.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Warning" -1 7 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 2 2 2 2 2 1 1 1 3 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 3 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "currentdir(\"C:\\\\ARCH\"):" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA 7 "" 1 "" {TEXT -1 50 "Warning, the name changecoords has been redefined\n " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "with(LinearAlgebra): wi th(Statistics): with(stats); with(stats[statplots]); with(RandomTools) :" }}{PARA 7 "" 1 "" {TEXT -1 40 "Warning, the name Rank has been rebo und\n" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7*%&anovaG%)describeG%$fitG%+ importdataG%'randomG%*statevalfG%*statplotsG%*transformG" }}{PARA 11 " " 1 "" {XPPMATH 20 "6#7.%(boxplotG%*histogramG%,scatterplotG%'xscaleG% 'xshiftG%+xyexchangeG%+xzexchangeG%'yscaleG%'yshiftG%+yzexchangeG%'zsc aleG%'zshiftG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "#with(Vect orCalculus):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 383 "sqIdentity Matrix := proc(n)\n description \"Returns an n-squared identi ty matrix\"; \n local M,i,j:\n M := Matrix(n,n):\n \+ for i from 1 to n do\n for j from 1 to n do\n \+ if evalb(i=j) then \n M[i, j] := 1:\n end if:\n end do:\n \+ end do:\n M\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "#Powells Method, \"Numerical Recipes in C\" p415" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2123 "Powells := proc(f,X,n)\n \+ description \"Find minimum of an n dimensional function f which \+ accepts an argument list of size n.\";\n local X_hat,fArgs,i,k, j,P,U,gamma_hat, gamma_hat_k:\n P := Array(0..n):\n X_ha t := Array(0..n):\n X_hat := X:\n # create n*n identity \+ matrix, init directions to basis vectors\n U := sqIdentityMatri x(n):\n for k from 0 to n do P[k] := Vector[row](1..n): end do: \n i := 0:\n while i < 10 do\n # starting pos \+ p0\n P[0] := X_hat:\n for k from 1 to n do \n \+ # construct f argument list\n fArgs := []:\n \+ for j from 1 to n do\n fArgs := [op(fAr gs), P[k-1][j]+gamma_hat*U[k,j]]:\n end do:\n \+ # move P[k-1] to the min of f(),along direction, U[k] \n \+ minimize(f(fArgs), gamma_hat, location):\n gamm a_hat_k := subs(%[2][1][1],gamma_hat):\n for j from 1 to n do\n P[k][j] := P[k-1][j] + gamma_hat_k*U[k,j]: \+ \n end do:\n end do:\n # set U[k] to U[k+1]\n for k from 1 to n-1 do \n for j from 1 to n do \n U[k,j] := U[k+1,j]:\n en d do:\n end do:\n # set U[N] to P[N]-P[0]\n \+ for j from 1 to n do \n U[n,j] := P[n][j] - P[0][j]: \n end do:\n i := i + 1:\n # construct f argument list\n fArgs := []:\n for j from 1 to n \+ do\n fArgs := [op(fArgs), P[0][j]+gamma_hat*U[n,j]]:\n \+ end do:\n # move P[n] to min, gamma_hat_k in direct ion u[n] and set to P[0] (x_hat)\n minimize(f(fArgs), gamma_ hat, location):\n gamma_hat_k := subs(%[2][1][1],gamma_hat): \n for j from 1 to n do\n X_hat[j] := P[0][j]+ gamma_hat_k*U[n,j]:\n end do:\n print(evalf(X_hat) );\n end do:\n fArgs := []:\n for j from 1 to n d o\n fArgs := [op(fArgs), X_hat[j]]:\n end do:\n \+ print(\"Minimum\", f(fArgs));\n X_hat\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 154 "plot_series := proc(A,B,a)\n \+ description \"Plot a time series\";\n pointplot([seq([t-A ,a[t]], t=A+1..B)], style=line, color=black)\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 264 "plot_series2 := proc(A,B,a,b)\n \+ description \"Plot two time series on the one graph\";\n \+ display(\{pointplot([seq([t-A,a[t]], t=A..B)], style=line, color=bl ack),pointplot([seq([t-A,b[t]], t=A..B)], style=line, color=COLOR(RGB, 0.7,0.7,0.7))\})\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 609 "AR := proc(T,P,vAlpha, Errors)\n description \"An autoregr essive series of length T\";\n local data,t,p,ar,e:\n da ta := array(1..T):\n ar := array(1..T):\n e := array(1.. T):\n for t from 1 to T do \n e[t] := stats[random, \+ normald[0,1]]():\n data[t] := 0:\n if t>P then \+ \n data[t] := vAlpha[1]*ar[t-1]: \n for p fr om 2 to P do \n data[t] := data[t] + vAlpha[p]*ar[t-p ]: \n end do: \n end if: \n \+ ar[t] := data[t]+e[t]:\n end do:\n ar \nend \+ proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 614 "MA := proc(T,P,vA lpha, Errors)\n description \"A Moving Average series of length T\";\n local data,t,p,ma,e:\n data := array(1..T):\n \+ ma := array(1..T):\n e := array(1..T):\n for t from 1 to T do \n e[t] := stats[random, normald[0,1]]():\n \+ data[t] := e[t]:\n if t>P then \n data[ t] := data[t] + vAlpha[1]*e[t-1]: \n for p from 2 to P do \n data[t] := data[t] + vAlpha[p]*e[t-p]: \n \+ end do: \n end if: \n ma[t ] := data[t]:\n end do:\n ma \nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 563 "ARCH := proc(T,P,vAlpha)\n \+ description \"Generates T Zero Mean ARCH(P) data\";\n local dat a,e,t,p,var:\n e := array(1..T):\n var := array(1..T):\n data := array(1..T):\n for t from 1 to T do\n \+ e[t] := stats[random, normald[0,1]]():\n var[t] := vAlpha [1]:\n if t > P then \n for p from 2 to P do \n var[t] := var[t] + vAlpha[p]*(data[t-(p-1)]^2):\n \+ end do:\n end if:\n data[t] := sqr t(var[t])*e[t]:\n end do:\n data\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 803 "GARCH := proc(T,P,Q,vAlpha, vBeta)\n description \"Generates T Zero Mean ARCH(P) data\";\n \+ local data,e,t,p,q,var:\n e := array(1..T):\n var := array(1..T):\n data := array(1..T):\n for t from 1 t o T do\n e[t] := stats[random, normald[0,1]]():\n \+ var[t] := 0:\n if P > 0 then var[t] := var[t] + vAlpha[1] : end if:\n if t > P then \n for p from 2 to \+ P do\n var[t] := var[t] + vAlpha[p]*(data[t-(p-1)]^2) :\n end do:\n end if:\n if t > Q t hen \n for q from 1 to Q do\n var[t] := var[t] + vBeta[q]*(var[t-q]^2):\n end do:\n \+ end if:\n data[t] := sqrt(var[t])*e[t]:\n end do:\n var \nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 168 "GenerateErrors := proc(T)\n local e,t:\n e := array(1 ..T):\n for t from 1 to T do\n e[t] := stats[random, norma ld[0,1]]():\n end do:\n e\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 497 "ARCH2 := proc(T,P,vAlpha,Errors)\n descr iption \"Generates T Zero Mean ARCH(P) data\";\n local data,t,p ,var:\n var := array(1..T):\n data := array(1..T):\n \+ for t from 1 to T do\n var[t] := vAlpha[1]:\n \+ if t > P then \n for p from 2 to P do\n \+ var[t] := var[t] + vAlpha[p]*(data[t-(p-1)]^2):\n end do:\n end if:\n data[t] := sqrt(var[t])*Errors[ t]:\n end do:\n data\nend proc:" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 53 "# p > 1 (so if p = 2, then really estimating A RCH(1) " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1433 "ARCH_estimate \+ := proc(T0,T,y,p)\n description \"Estimates a zero mena ARCH(p-1) time series for observed data y of length T. Returns an ordered list \+ of the parameter estimates for alpha\":\n local dL, alphaInterval ,h,L,z,delta,alpha,i,solution,estimates:\n alpha := 'alpha':\n \+ # create paramater vector z=1,y[t-1]^2,y[t-2]^2,... delta=alpha[1], alpha[2],... \n z := Vector(p):\n delta := Vector(p):\n z[1] := 1:\n delta[1] := alpha[1]:\n for i from 2 to p do\n z[i] := y[t-(i-1)]^2:\n delta[i] := alpha[i]: \n end do:\n print(Transpose(z)):\n print(delta):\n \+ print(Multiply(Transpose(z), delta)): \n # arch disturbance val ues h = [z'][delta]\n h := (alpha) -> Multiply(Transpose(z),delta );\n # log likelihood\n L := (alpha) -> (T/2)*log(2*Pi) + (T /2)*Sum(log(h(alpha)),t=p..T) + (T/2)*Sum((y[t]^2)/h(alpha),t=p..T);\n dL := \{\}:\n alphaInterval := \{\}:\n # calc derivati ves of L\n for i from 1 to p do\n dL := \{op(dL),value(D iff(L(alpha),alpha[i]))=0\}:\n alphaInterval := \{op(alphaInt erval),alpha[i]=0..1\}:\n end do: \n # print(dL):\n \+ print(alphaInterval):\n solution := value(fsolve(dL,alphaInterv al)):\n estimates := []:\n # get alpha parameter estimates\n for i from 1 to p do\n estimates := [op(estimates),sub s(solution, alpha[i])];\n end do:\n estimates\nend proc:" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 254 " AIC := proc(T0,T,y,alpha) \n description \"Akaike's IC for zero mean ARCH(p) model using ML parameter estimates alpha and observed data y. p is size of alpha\": \n local p:\n p := op(1,alpha):\n value(evalf(-L(T0,T,y ,alpha)+2*p/T))\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 274 "AICc := proc(T0,T,y,alpha)\n description \"Akaike's Correcte d IC for zero mean ARCH(p) model using ML parameter estimates alpha an d observed data y. p is size of alpha\":\n local p:\n p := o p(1,alpha):\n value(evalf(-L(T0,T,y,alpha)+2*(p+1)/(T-p-2)))\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 270 "BIC := proc(T0,T ,y,alpha)\n description \"Schwart's Bayesian IC for zero mean ARC H(p) model using ML parameter estimates alpha and observed data y. p i s size of alpha\":\n local p:\n p := op(1,alpha):\n val ue(evalf(-L(T0,T,y,alpha)+(p*log(T))/T))\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 276 "HQ := proc(T0,T,y,alpha)\n descript ion \"Hannan and Quinn's IC for zero mean ARCH(p) model using ML param eter estimates alpha and observed data y. p is size of alpha\":\n \+ local p:\n p := op(1,alpha):\n value(evalf(-L(T,T0,y,alpha) +(2*p*log(log(T)))/T))\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 526 "likelihood := proc(T0,T,y,alpha_hat) \n local i ,p,z,alpha,Alpha,h,L,y_var:\n p := op(1,alpha_hat):\n z := V ector(p):\n z[1] := 1:\n alpha := 'alpha':\n Alpha := [ alpha[1]]:\n for i from 2 to p do\n z[i] := y[t-(i-1)]^2 :\n Alpha := [op(Alpha),alpha[i]]:\n end do:\n h := (Alpha) -> Multiply(Transpose(z),Vector[column](Alpha));\n L := \+ (Alpha) -> Product( (1/(sqrt(2*Pi*h(Alpha)))) * exp( (-y[t]^2) / (2*h( Alpha)) ), t=T0+1..T+T0);\n value(evalf(L(alpha_hat)))\n end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 662 "L := proc(T0,T,y,alpha_h at) \n local i,p,z,alpha,Alpha,h,L,y_var:\n p := op(1,alpha_ hat):\n z := Vector(p):\n z[1] := 1:\n alpha := 'alpha' :\n Alpha := [alpha[1]]:\n for i from 2 to p do\n z [i] := y[t-(i-1)]^2:\n Alpha := [op(Alpha),alpha[i]]:\n \+ end do:\n h := (Alpha) -> Multiply(Transpose(z),Vector[column](Al pha));\n # first h(alpha) for t=1..p-1 should be unconditional va riance\n #y_var := stats[describe,variance](convert(y,list)):\n \+ L := (Alpha) -> -(T/2)*log(2*Pi) - (1/2)*Sum(log(h(Alpha)),t=T0+1. .T+T0) - (1/2)*Sum((y[t]^2)/h(Alpha),t=T0+1..T+T0);\n value(evalf(L(a lpha_hat)))\n end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 624 "Fisher := proc(T0,T,y,alpha_hat) \n local i,p,z,alpha,Alpha, h,L,y_var,Info:\n p := op(1,alpha_hat):\n z := Vector(p):\n \+ z[1] := 1:\n alpha := 'alpha':\n Alpha := [alpha[1]]:\n for i from 2 to p do\n z[i] := y[t-(i-1)]^2:\n \+ Alpha := [op(Alpha),alpha[i]]:\n end do:\n h := (Alpha) -> \+ Multiply(Transpose(z),Vector[column](Alpha));\n # first h(alpha) \+ for t=1..p-1 should be unconditional variance\n Info := (Alpha) - > evalf(value(1/(2*T)*Sum(Multiply(Matrix(z),Matrix(Transpose(z))/h(Al pha)^2), t=T0..T+T0))):\n value(evalf(Determinant(Info(alpha_hat) )))\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1446 "ARCH_sc ore := proc(T0,T,y,alpha_0,N) \n description \"Estimate ARCH para meters with alpha_0 initial values\":\n local i,n,p,z,alpha,newal pha, Alpha,h,l,dl,Info,InfoInverse,grad,dalpha,y_var: \n p := nops(alpha_0):\n z := Vector(p):\n z[1] := 1:\n alpha \+ := 'alpha':\n Alpha := [alpha[1]]:\n for i from 2 to p do\n \+ z[i] := y[t-(i-1)]^2:\n Alpha := [op(Alpha),alpha[i] ]:\n end do:\n h := (Alpha) -> Multiply(Transpose(z),Vector[ column](Alpha));\n # kickstart\n y_var := stats[describe,var iance](convert(y,list)):\n l := (Alpha) -> -(1/2)*log(2*Pi) - (1/ 2)*log(h(Alpha)) - (1/2)*(y[t]^2)/h(Alpha);\n\n dl := []:\n \+ for i from 1 to p do\n dl := [op(dl),value(Diff(l(Alpha),alph a[i]))]:\n end do: \n alpha := Vector(alpha_0):\n n := \+ 0:\n while evalb(n " 0 "" {MPLTEXT 1 0 562 "VarCov := \+ proc(T,y,alpha_hat) \n local z,Alpha,alpha,i,p,h,LL,Info,mml, pri or_alpha, prior_p:\n p := op(1,alpha_hat):\n z := Vector(p): \n z[1] := 1:\n alpha := 'alpha':\n Alpha := [alpha[1]] :\n for i from 2 to p do\n z[i] := y[t-(i-1)]^2:\n \+ Alpha := [op(Alpha),alpha[i]]:\n end do:\n h := (Alpha) \+ -> Multiply(Transpose(z),Vector[column](Alpha));\n Info := (Alpha ) -> evalf(value(1/(2*(T-(p-1)))*Sum(Multiply(Matrix(z),Matrix(Transpo se(z))/h(Alpha)^2), t=p..T))):\n MatrixInverse(Info(alpha_hat))\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 223 "EstimatesZero := proc(A)\n local i,count:\n count := 0:\n for i from 1 to op(1,A) do\n if(evalb(A[i] = 0)) then\n c ount := count + 1:\n end if:\n end do:\n count\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1528 "MML87 := proc(T 0,T,y,alpha_hat,prior_num,small_sample_num)\n local z,Alpha,alpha ,i,p,h,LL,Info,mml, prior_alpha, prior_p, lattice, D:\n p := op(1 ,alpha_hat):\n prior_p := 1:\n \n\n if(evalb(prior_num = 1)) then\n if (evalb(p > 1)) then \n prior_alph a := 1/(.99/(p-1))^(p-1); \n else\n prior_alpha := \+ 1:\n end if:\n elif(evalb(prior_num = 2)) then\n \+ if (evalb(p > 1)) then \n prior_alpha := ((1-add(alpha_ hat[i], i=2..p))/max(0.0000000001,alpha_hat[1]) )^(1/2);\n el se\n prior_alpha := 1:\n end if:\n elif(eval b(prior_num = 3)) then\n if (evalb(p > 2)) then \n \+ prior_alpha := (p-2)!:\n else \n prior_alpha := 1: \n end if:\n else\n print( \"FUCK!\");\n end if:\n\n if (evalb(p > 3)) then\n \+ lattice := 1/2*Pi*exp(1);\n elif(evalb(p = 3)) then \n \+ lattice := 5/(36*sqrt(3)):\n elif(evalb(p = 2)) then\n \+ lattice := 19/(192*2^(1/3)):\n elif(evalb(p = 1)) then\n \+ lattice := 1/12:\n end if:\n\n if(evalb(small_sample _num = 0)) then\n mml := -log(prior_alpha/(sqrt(Fisher(T0,T,y,a lpha_hat)/12)))-L(T0,T,y,alpha_hat)+log(p/12)+p/2:\n elif(evalb(s mall_sample_num = 1)) then\n mml := - (1/2)*log(1+(Fisher(T0,T, y,alpha_hat)*lattice)/prior_alpha^2)- L(T0,T,y,alpha_hat) + log(p/12)+ p/2:\n else\n print(\"Fuck\");\n end if:\n\n value (evalf(mml))\nend proc: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 233 "ExportArray := proc(filename, data) \n description \+ \"write data to currentdir\\data\\filename in vector column format\": \n ExportVector(cat(currentdir(),cat(\"\\\\\",filename)), convert(data,Vector[column])):\nend proc:" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 213 "ImportArray := proc(filename) \n desc ription \"read data from currentdir\\data\\filename and return as a li st\": \n convert(ImportVector(cat(currentdir(),cat(\"\\\\ \",filename))),'list')\n end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 380 "# Note: overwrites existing files \ngenerateData := \+ proc(N,T,p,alpha) \n description \"generate N sets of s ize T ARCH(p) data with alpha parameters and write to files archp-T-i. dat\":\n local i, data:\n for i from 1 t o N do \n ExportData(cat(\"arch\",p,\"-\",T,\"-\",i ,\".dat\"), ARCH(T,p,alpha)):\n end do:\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2772 "order_choice := proc(crite rion,T0,T,y,estimates,pmax)\n #description \"Return the mode l order chosen by the criterion from the paramter estimates for pmax m odels\":\n local i,current,best,bestp,mmldRet:\n\n \+ if(evalb(criterion = \"MML87a\")) then\n best := MML87(T 0,T,y,estimates[1],1,0):\n elif(evalb(criterion = \"MML87a_s ml\")) then\n best := MML87(T0,T,y,estimates[1],1,1):\n \+ elif(evalb(criterion = \"MML87b\")) then\n best \+ := MML87(T0,T,y,estimates[1],2,0):\n elif(evalb(criterion = \+ \"MML87b_sml\")) then\n best := MML87(T0,T,y,estimates[1] ,2,1):\n elif(evalb(criterion = \"MML87c\")) then\n \+ best := MML87(T0,T,y,estimates[1],3,0):\n elif(evalb(cr iterion = \"MML87c_sml\")) then\n best := MML87(T0,T,y,es timates[1],3,1):\n elif(evalb(criterion = \"BIC\")) then\n \+ best := BIC(T0,T,y,estimates[1]):\n elif(evalb(c riterion = \"AIC\")) then\n best := AIC(T0,T,y,estimates[ 1]):\n elif(evalb(criterion = \"AICc\")) then\n \+ best := AICc(T0,T,y,estimates[1]):\n elif(evalb(criterion = \"HQ\")) then\n best := HQ(T0,T,y,estimates[1]):\n \+ end if:\n bestp := 1:\n\n for i from 2 to pma x do\n if(evalb(criterion = \"MML87a\")) then\n \+ current := MML87(T0,T,y,estimates[i],1,0):\n elif(evalb(c riterion = \"MML87b\")) then\n current := MML87(T0,T,y,es timates[i],2,0):\n elif(evalb(criterion = \"MML87c\")) then \n current := MML87(T0,T,y,estimates[i],3,0):\n \+ elif(evalb(criterion = \"MML87a_sml\")) then\n current \+ := MML87(T0,T,y,estimates[i],1,1):\n elif(evalb(criterion = \"MML87b_sml\")) then\n current := MML87(T0,T,y,estimate s[i],2,1):\n elif(evalb(criterion = \"MML87c_sml\")) then\n current := MML87(T0,T,y,estimates[i],3,1):\n \+ elif(evalb(criterion = \"BIC\")) then\n current := BIC(T0 ,T,y,estimates[i]):\n elif(evalb(criterion = \"AIC\")) then \n current := AIC(T0,T,y,estimates[i]):\n elif (evalb(criterion = \"AICc\")) then\n current := AICc(T0,T ,y,estimates[i]):\n elif(evalb(criterion = \"HQ\")) then\n \+ current := HQ(T0,T,y,estimates[i]):\n end if:\n \n if(evalb(current < best)) then\n be st := current:\n bestp := i:\n end if: \n end do: \n\n # return criterion u sed, model order (num parameters inc. alpha[0]?) selected, with -LogLi kelihood for the model selected\n [criterion, op(1,estimates[be stp]),-L(T0,T,y,estimates[bestp])] end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 893 "MMLD_order_choice := proc(T0,T,y,estimates,pmax ,SampleSize)\n #description \"Return the model order chosen by the criterion from the paramter estimates for pmax models\":\n \+ local i,current,best,bestp,mmldRet:\n\n best := MMLD( T0,T,y,estimates[1],SampleSize);\n bestp := 1:\n\n \+ for i from 2 to pmax do\n #print(\"MMLD \", i);\n\n \+ current := MMLD(T0,T,y,estimates[i],SampleSize):\n \+ if(evalb(current[2] < best[2])) then\n best := cu rrent:\n bestp := i:\n end if: \+ \n end do: \n\n print(best[1]);\n # \+ return criterion used, model order (num parameters inc. alpha[0]?) sel ected, with -LogLikelihood for the model selected\n # and param eter estimates\n [\"MMLD\", nops(best[1]),-L(T0,T,y,convert(bes t[1],Vector)),best[1]]\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1053 "getStats := proc(criterion,T0,T,y,errors,p,pmax,est imates) \n description \"Returns stats p,avg -LogLikeli hood \n for model criterion\": \n local ret, stats,mspe, \+ y_hat:\n stats := [0,0,0,0,0]:\n \n if(evalb(cr iterion = \"MMLD\")) then\n ret := MMLD_order_choice(T0,T,y, estimates,pmax,200):\n else \n ret := order_choice(c riterion,T0,T,y,estimates,pmax);\n end if:\n\n if(eval b(ret[2]p)) then \n stats[3] := stats[3]+1:\n else\n print(\"FCUK\"):\n end if:\n\n \+ stats[4] := stats[4] + ret[3]; # get log likelihood \n\n \+ if(evalb(criterion = \"MMLD\")) then\n stats[5] := MSPE(T0, T,y,errors,nops(ret[4]),ret[4]); # CHECK nops(ret[1]) is corre ct\n else\n stats[5] := MSPE(T0,T,y,errors,ret[2],e stimates[ret[2]]):\n end if:\n stats \nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 364 "MSPE := proc(T0,T,y,e,p,alp ha_hat)\n local t,z,i,Alpha,h,mspe,alpha,y_hat;\n mspe := 0: \n for t from T0 to T+T0 do\n h := alpha_hat[1]:\n \+ for i from 2 to p do \n h := h + alpha_hat[i]*y[ t-(i-1)]^2:\n end do:\n mspe := mspe + (y[t] - e [t]*sqrt(h))^2; \n end do:\n mspe/T\n end proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 151 " ExportAlphaEstimates := proc(filename,estimates)\n ExportVector(cat(currentdir(),cat(\" \\\\\",filename)),convert(estimates,Vector[column])):\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 549 "ImportAlphaEstimates := pro c(filename) \n local i,j,k,filedata,alpha,alpha_estimates:\n \+ k := 1:\n i := 1: \n filedata := ImportVector(cat(currentdi r(),cat(\"\\\\\", filename)));\n alpha_estimates := []:\n wh ile evalb(k " 0 "" {MPLTEXT 1 0 402 "RandomAlpha := proc(p)\n local i, al pha,upper:\n alpha := Array(1..p):\n alpha[1] := stats[random, uniform[0, 1]]():\n if(evalb(p > 1)) then \n alpha[2] := stats[random, uniform[0, 1]]();\n for i from 3 to p do\n \+ upper := 1-add(alpha[j],j=2..(i-1)):\n alpha[i] := stats[random, uniform[0,upper]]():\n end do:\n end if :\n alpha\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 516 "GenerateTrials2 := proc(N,T0,T,p,alpha) \n description \"Gene rate N trials of T ARCH(p) data with parameters alpha\";\n local i ,n,y,errors:\n #start trials\n for n from 1 to N do\n \+ #print('Trial', n);\n\n # generate Vt data\n errors := GenerateErrors(T+T0):\n ExportArray(cat(\"arch\",p,\"-\",T,\" -\",n,\"-errors.dat\"),errors):\n\n # generate true data\n \+ y := ARCH(T+T0,p,alpha): \n ExportArray(cat(\"arch \",p,\"-\",T,\"-\",n,\".dat\"),y):\n end do: \nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 612 "GenerateTrials := proc(N,T0 ,T,p) \n description \"Generate N trials of T ARCH(p) data with ra ndom parameters alpha, kickstart with T0 values\";\n local i,n,y,a lpha, errors:\n #start trials\n for n from 1 to N do\n \+ #print('Trial', n);\n\n # get 'true' parameters\n alp ha := RandomAlpha(p);\n\n # generate Vt data\n errors \+ := GenerateErrors(T+T0):\n ExportArray(cat(\"arch\",p,\"-\",T, \"-\",n,\"-errors.dat\"),errors):\n\n # generate true data\n \+ y := ARCH2(T+T0,p,alpha,errors):\n ExportArray(cat(\"ar ch\",p,\"-\",T,\"-\",n,\".dat\"),y):\n end do:\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 672 "EstimateTrials := proc(N,T0 ,T,p,pmax) \n local i,n,alpha_estimates,alpha0,y:\n for n from 1 to N do\n #printf(\"Estimating %d...\", n);\n# print( cat(\"Estimating \",n,\"...\"));\n # import y\n y := Imp ortArray(cat(\"arch\",p,\"-\",T,\"-\",n,\".dat\")):\n # build i nitial alpha estimates, then estimate\n alpha_estimates := []: \n alpha0 := []:\n for i from 1 to pmax do\n \+ alpha0 := [op(alpha0), 0.5]:\n alpha_estimates := [op(a lpha_estimates), ARCH_score(T0,T,y,alpha0,5)]; \n end do: \n ExportAlphaEstimates(cat(\"arch\",p,\"-\",T,\"-\",n,\"-estim ates.dat\"),alpha_estimates):\n end do:\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 678 "EstimateTrials2 := proc(Na,Nb,T0,T ,p,pmax) \n local i,n,alpha_estimates,alpha0,y:\n for n from N a to Nb do\n #printf(\"Estimating %d...\", n);\n# print(c at(\"Estimating \",n,\"...\"));\n # import y\n y := Impo rtArray(cat(\"arch\",p,\"-\",T,\"-\",n,\".dat\")):\n # build in itial alpha estimates, then estimate\n alpha_estimates := []:\n alpha0 := []:\n for i from 1 to pmax do\n \+ alpha0 := [op(alpha0), 0.5]:\n alpha_estimates := [op(alp ha_estimates), ARCH_score(T0,T,y,alpha0,5)]; \n end do:\n ExportAlphaEstimates(cat(\"arch\",p,\"-\",T,\"-\",n,\"-estimat es.dat\"),alpha_estimates):\n end do:\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 4226 "SelectModels := proc(N,T0,T,p,pmax) \n \+ local y_filename, BIC_stats,AIC_stats,AICc_stats,HQ_stats,MMLD_sta ts, n,y,alpha_estimates,errors, justMML, MML87a_stats,MML87b_stats, MM L87c_stats, \n MML87a_sml_stats,MML87b_sml_stats, MML87c_sm l_stats;\n #init stats\n BIC_stats := [0,0,0,0,0]:\n AIC_sta ts := [0,0,0,0,0]:\n AICc_stats := [0,0,0,0,0]:\n HQ_stats := \+ [0,0,0,0,0]:\n MML87a_stats := [0,0,0,0,0]:\n MML87b_stats := \+ [0,0,0,0,0]:\n MML87c_stats := [0,0,0,0,0]:\n MML87a_sml_stats := [0,0,0,0,0]:\n MML87b_sml_stats := [0,0,0,0,0]:\n MML87c_s ml_stats := [0,0,0,0,0]:\n\n MMLD_stats := [0,0,0,0,0]:\n\n ju stMML = 1: #ie. 1=true=calc only mml \n\n #start trials\n prin t(cat(\"arch\",p,\"-\",T));\n for n from 1 to N do\n # impo rt y\n y_filename := cat(\"arch\",p,\"-\",T,\"-\",n,\".dat\"): \n y := ImportArray(y_filename):\n\n errors := ImportArr ay(cat(\"arch\",p,\"-\",T,\"-\",n,\"-errors.dat\"));\n alpha_es timates := ImportAlphaEstimates(cat(\"arch\",p,\"-\",T,\"-\",n,\"-esti mates.dat\")):\n# print(y_filename);\n # get stats for thi s trial\n\n if(evalb(justMML=1)) then\n BIC_stats := BIC _stats + getStats(\"BIC\",T0,T,y,errors,p,pmax,alpha_estimates):\n \+ AIC_stats := AIC_stats + getStats(\"AIC\",T0,T,y,errors,p,pmax, alpha_estimates):\n AICc_stats := AICc_stats + getStats(\"AI Cc\",T0,T,y,errors,p,pmax,alpha_estimates):\n HQ_stats := HQ _stats + getStats(\"HQ\",T0,T,y,errors,p,pmax,alpha_estimates):\n \+ MML87a_stats := MML87a_stats + getStats(\"MML87a\",T0,T,y,errors,p,p max,alpha_estimates):\n MML87b_stats := MML87b_stats + getStats( \"MML87b\",T0,T,y,errors,p,pmax,alpha_estimates):\n # MML87c_stat s := MML87c_stats + getStats(\"MML87c\",T0,T,y,errors,p,pmax,alpha_est imates):\n MML87a_sml_stats := MML87a_sml_stats + getStats(\"MML 87a_sml\",T0,T,y,errors,p,pmax,alpha_estimates):\n MML87b_sml_st ats := MML87b_sml_stats + getStats(\"MML87b_sml\",T0,T,y,errors,p,pmax ,alpha_estimates):\n # MML87c_sml_stats := MML87c_sml_stats + get Stats(\"MML87c_sml\",T0,T,y,errors,p,pmax,alpha_estimates):\n en d if:\n\n MMLD_stats := MMLD_stats + getStats(\"MMLD\",T0,T,y,er rors,p,pmax,alpha_estimates):\n\n # Update stats as we go along a nd write to file ( write n say we can change the stats later\n w ritedata[append](\"stats_mmld.dat\",[AIC_stats,AICc_stats,BIC_stats,HQ _stats,MML87a_stats,MML87b_stats,MML87a_sml_stats,MML87b_sml_stats,MML D_stats]);\n writedata(cat(\"arch\",p,\"-\",T,\"-\",n,\"-stats_m mld.dat\"),[AIC_stats,AICc_stats,BIC_stats,HQ_stats,MML87a_stats,MML87 b_stats,MML87a_sml_stats,MML87b_sml_stats,MMLD_stats]);\n\nend do:\n\n \n # average -L\nif(evalb(justMML=1)) then\n BIC_stats[4] := B IC_stats[4]/N:\n AIC_stats[4] := AIC_stats[4]/N:\n AICc_stats[ 4] := AICc_stats[4]/N:\n HQ_stats[4] := HQ_stats[4]/N:\n MML87 a_stats[4] := MML87a_stats[4]/N:\n MML87b_stats[4] := MML87b_stats [4]/N:\n # MML87c_stats[4] := MML87c_stats[4]/N:\n MML87a_sml_s tats[4] := MML87a_sml_stats[4]/N:\n MML87b_sml_stats[4] := MML87b_ sml_stats[4]/N:\n # MML87c_sml_stats[4] := MML87c_sml_stats[4]/N:\n end if:\n MMLD_stats[4] := MMLD_stats[4]/N:\n\n # average MSPE \nif(evalb(justMML=1)) then\n BIC_stats[5] := BIC_stats[5]/N:\n \+ AIC_stats[5] := AIC_stats[5]/N:\n AICc_stats[5] := AICc_stats[5] /N:\n HQ_stats[5] := HQ_stats[5]/N:\n MML87a_stats[5] := MML87 a_stats[5]/N:\n MML87b_stats[5] := MML87b_stats[5]/N:\n # MML87 c_stats[5] := MML87c_stats[5]/N:\n MML87a_sml_stats[5] := MML87a_s ml_stats[5]/N:\n MML87b_sml_stats[5] := MML87b_sml_stats[5]/N:\n \+ # MML87c_sml_stats[5] := MML87c_sml_stats[5]/N:\nend if:\n\n MMLD_st ats[5] := MMLD_stats[5]/N:\n\n # output stats\nif(evalb(justMML=1) ) then\n print('bic', BIC_stats);\n print('aic', AIC_stats);\n print('aicc', AICc_stats);\n print('hq', HQ_stats);\n prin t('mml87a', MML87a_stats);\n print('mml87b', MML87b_stats);\n pr int('mml87a_sml', MML87a_sml_stats);\n print('mml87b_sml', MML87b_s ml_stats);\nend if:\n# print('mmld', MMLD_stats);\n writedata(\"r esults_auto_mmld.dat\",[AIC_stats,AICc_stats,BIC_stats,HQ_stats,MML87a _stats,MML87b_stats,MML87a_sml_stats,MML87b_sml_stats,MMLD_stats]);\ne nd proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 572 "MetropolisHast ings := proc(T,T0,f,q,Q,x0)\n description \"sample T points from \+ distribution f using proposal dist q (and Q)\";\n local t,i,X,U,Y ,AcceptProb:\n AcceptProb := (X,Y) -> min(1,(f(Y)*q(X,Y))/(f(X)*q (Y,X)));\n t := 1:\n X := Array(1..T):\n X[1] := x0:\n \+ for i from 1 to T-1 do\n U := stats[random, uniform[0,1] ]():\n Y := Q(X[t]):\n if(evalb(U <= AcceptProb(X[t] ,Y))) then\n X[t+1] := Y:\n else\n X[ t+1] := X[t]:\n end if:\n t := t+1:\n end do:\n X\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 136 "Burn In := proc(t0,T,X) \n local data:\n data := Array(1. .T-t0):\n Copy(T-t0,X,t0,data);\n data\nend proc:" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 425 "graphMCMC := proc(T,x) \n local mean, var, sd, pUB,pLB,pMean, pData:\n pData := pointpl ot([seq([t,x[t]], t=1..T)], style=line, color=black);\n mean := st ats[describe,mean](convert(x,list));\n var := stats[describe,vari ance](convert(x,list));\n sd := sqrt(var):\n display(\{pData,p lot(mean,1..T), plot(mean+2*sd,1..T,linestyle=DASH), plot(mean-2*sd,1. .T,linestyle=DASH)\}, axes=boxed, title=`MCMC`);\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 183 "SortSample := proc(P,L,N) \+ \n# Sorts the sample points P according to L\n# returns as a list (not a array)\nconvert(InsertionSort(convert(L,array),convertSampleToVec(P ),N),list)\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 142 "c onvertSampleToVec := proc(S) \nlocal vecS,i:\nvecS := array(1..nops(S) ):\nfor i from 1 to nops(S) do\n vecS[i] := S[i]:\nend do:\nvecS\ne nd proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 477 "InsertionSort \+ := proc(A,B,N)\n#insertion sort, user defined to keep pointers to wher e each element is positioned\n# sorts array A but also keeps B in the \+ same order as A \nlocal i,j,AA,BB,tmpA,tmpB:\nAA := A:\nBB := B:\nfor \+ i from 1 to N do\n tmpA := AA[i]:\n tmpB := BB[i]:\n j := i: \n while (j > 1 and AA[j-1] < tmpA) do\n AA[j] := AA[j-1]: \n BB[j] := BB[j-1]:\n j := j - 1:\n end do: \n AA[j] := tmpA: \n BB[j] := tmpB:\nend do:\nBB\nend proc:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 472 "SampleExp := proc(N,p,vB,vL amda) \n #Generate N random p-dimensional points from p-dimensional exponential dist centred at B vector with lamda vector\n local pts ,i,j,B,lamda,EXP,E:\n EXP := Distribution(PDF = (x -> piecewise(x<0 ,0,lamda*exp((B-x)*lamda)))):\n E := RandomVariable(EXP):\n \n p ts := array(1..p):\n for i from 1 to p do\n lamda := vLamda[i] : B := vB[i]: pts[i] := Sample(E,N):\n end do:\n [seq( [seq(pts[ i][j], i=1..p)], j=1..N)]\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 362 "pdfExp := proc(vB,vLamda,p,pt)\n local i,prod, E,EXP,B,lamda:\n EXP := Distribution(PDF = (x -> piecewise(x<0,0 ,lamda*exp((B-x)*lamda)))):\n E := RandomVariable(EXP):\n \+ prod := 1:\n for i from 1 to p do\n B := vB[i]: lamda \+ := vLamda[i]:\n prod := prod * PDF(E,pt[i],numeric):\n \+ end do:\n evalf(prod)\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 594 "checkEnvelope := proc(f,g,A,B,step) \n des cription \"Returns false if posterior f([i,j]) exceeds envelope g([i,j ])\":\n local i,j,df,ret,exceed:\n exceed = false: \n ret := []:\n for i from A to B by step do \n \+ for j from A to B by step do\n if(evalb( g([i,j])-f([i,j]) <= 0)) then \n exceed = true:\n ret := [i,j,f([i,j])]:\n i \+ := B: j := B: \n break:\n end i f:\n end do:\n end do:\n ret\nend proc: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 3537 "MMLD := proc(T0,T,y,a lpha_hat,N)\n local i,j,U,p,alpha_ml,h,f,logf,posterior,envelop e_ml, pts, vB,vLamda,S,Q,maxS,ratio,\n nit,firstNumer,fir stDenom,secondNumer,secondDenom,secondLength,Sf,Ssorted,oneOnF,msglen, pointEstimate:\n\n print(\"Running MMLD\");\n\n p := op( 1,alpha_hat):\n if (evalb(p > 1)) then \n h := \+ 1/((.99/(p-1))^(p-1)):\n else\n h := 1: \n \+ end if:\n \n #Map f(y,alpha) -> libray likelihoo d function\n f := (y,alpha) -> likelihood(T0,T,y,convert(alpha, Vector)):\n logf := (y,alpha) -> L(T0,T,y,convert(alpha,Vector) ):\n\n posterior := (alpha) -> h(alpha)*logf(y,alpha):\n\n \+ envelope_ml := 0.9*posterior(alpha_hat):\n\n # generate poin ts from envelope dist\n vB := alpha_hat: \n vLamda := [] :\n for i from 1 to p do\n vLamda := [op(vLamda),1]: \n end do:\n\n # take bigger sample from envelope than N for random coding later\n pts := SampleExp(N*2,p,vB,vLamda):\n \n\n#Get S\n # Generate Sample S from posterior (using points f rom exponential dist as weights)\n ratio := []: # keep ratios f or later use\n maxS := posterior(pts[1])/(pdfExp(vB,vLamda,p,pt s[1])+envelope_ml):\n for i from 1 to N do\n ratio : = [op(ratio),posterior(pts[i])/(pdfExp(vB,vLamda,p,pts[i])+envelope_ml )]:\n# print(ratio[i]); \n if(evalb(ra tio[i]>maxS)) then \n maxS := ratio[i]:\n end if:\n end do:\n S := []:\n for i from 1 to N do\n \+ U := stats[random, uniform[0,1]]():\n if(evalb(U <= ratio[i]*(1/maxS))) then \n S := [op(S),pts[i]]:\n \+ end if: \n end do:\n\n#Get Q\n # Sort S in terms of l ikelihood\n Sf := []:\n for i from 1 to nops(S) do\n \+ Sf := [op(Sf), logf(y,S[i])]: \n end do:\n Ssorted := \+ SortSample(S,Sf,nops(Sf)):\n \n # Fitzgibbon (2004) Algorithm \+ 6.1\n # find subset Q of S S=Ssorted:\n i := 'i':\n \+ S := Ssorted:\n nit := evalf(log[2](exp(1))):\n firstNumer := 1/f(y,S[1]):\n firstDenom := add(1/f(y,S[i]),i=1..nops(S)): \n secondNumer := -firstNumer * logf(y,S[1]):\n secondDeno m := firstNumer:\n secondLength := secondNumer / secondDenom:\n \+ Q := []:\n Q := [op(Q),S[1]]:\n i := 2:\n whil e ( i < nops(S) and -logf(y,S[i]) <= secondLength + 1 ) do\n \+ oneOnF := 1/f(y,S[1]):\n firstNumer := firstNumer + oneOnF:\n \+ secondNumer := secondNumer - oneOnF * logf(y,S[i]):\n \+ secondDenom := secondDenom + oneOnF:\n secondLength := secondN umer / secondDenom:\n Q := [op(Q),S[i]]:\n i := i + 1: \n end do:\n\n# Calc MsgLen\n # MMLD Message Length\n \+ # Fitzgibbon (2004) Equation 6.9\n firstNumer := add(1/f(y,Q[i ]), i=1..nops(Q)):\n #firstDenom := add(1/f(y,S[i]), i=1..nops(S )): # USE firstDenom from when finding Q?\n secondNumer := add(- logf(y,Q[i])/f(y,Q[i]), i=1..nops(Q)): # could use f(y,Q[i]) then lo g(f(.)) instead of calculating f(.) twice\n secondDenom := first Numer:\n msglen := -log(firstNumer/firstDenom) + secondNumer/sec ondDenom + evalf(log(p)):\n \n# Get Point Estimate\n pointEstimat e := MMLDpointEstimate(T0,T,p,y,[op(pts[N..N*2])],secondLength+1,poste rior,vB,vLamda,envelope_ml);\n# pointEstimate := MMLDpointEstimat e(T0,T,p,y,pts,ratio,maxS,secondLength+1,posterior,vB,vLamda,envelope_ ml);\n\n [pointEstimate,msglen]\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 1521 "MMLDpointEstimate := proc(T0,T,p,y,pts,boundar y,posterior,vB,vLamda,envelope_ml) \n# Random Coding of Estimates \n# \+ Wallace (2005) Section 4.10.1\n\n local i,S,U,N,max_tries,tries,e stimate,n,s,ratio,maxS:\n # get first random element in S that is in the region R \n print(\"Estimating Points\");\n N := nop s(pts):\n \n # Generate Sample S from posterior (using points fr om exponential dist as weights)\n ratio := []: # keep ratios for \+ later use\n maxS := posterior(pts[1])/(pdfExp(vB,vLamda,p,pts[1]) +envelope_ml):\n for i from 1 to N do\n ratio := [op (ratio),posterior(pts[i],y)/(pdfExp(vB,vLamda,p,pts[i])+envelope_ml)]: \n if(evalb(ratio[i]>maxS)) then \n maxS := r atio[i]:\n end if:\n end do:\n\n# S := pts:\n S := []:\n for i from 1 to N do\n U := stats[ra ndom, uniform[0,1]]():\n if(evalb(U <= ratio[i]*(1/maxS))) t hen \n S := [op(S),pts[i]]:\n end if: \n e nd do:\n\n# S generated from posterior, need to weight by 1/likelihood to get from prior\n max_tries := N;\n estimate := 1:\n \+ for n from 1 to nops(S) do \n i := Generate(integer(range =1..nops(S))):\n s := S[i]/ likelihood(T0,T,y,convert(S[i] ,Vector));\n if(evalb(-L(T0,T,y,convert(s,Vector)) < bound ary)) then\n estimate := i:\n break: \+ \n end if:\n end do: \+ \n print(i): \n S[estimate]\nend proc:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "54 0 0" 567 }{VIEWOPTS 1 1 0 3 4 1802 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }