請教各位
如何建立一個規劃求解的程式ㄋ
我已經弄好幾天
但是一值抓不到頭緒
來這請各位幫忙
謝謝

Posted by vbqa at 痞客邦 PIXNET 留言(11) 引用(0) 人氣()


open trackbacks list Trackbacks (0)

留言列表 (11)

Post Comment
  • Shawn
  • 是哪一種方程式啊?
    二元一次? 聯立? 還是微分方程??
    你可以用矩陣方法求一般equation的解....
    for example:
    a1x+b1y=c1
    a2x+b2y=c2
    solution:
    x=[c1 b1
    c2 b2]/delta
    y=[a1 c1
    a2 c2]/delta
    where delta=[a1 b1
    a2 b2]

  • les
  • 首先謝謝Shawn您的回應
    我將程式的列出來與您討論
    希望您能試著幫忙我解決這個問題
    因為這關係到我能不能畢業
    min: Var(rp) = w1^2 * var1 + w2^2 * var2 + 2 * w1 * w2 * sig1 * sig2
    E(rp) = w1 * E(r1) + w2 * E(r2)
    其中 w1+w2 = 1
    var1 = sig1^2
    var2 = sig2^2
    其實變數與權重的部分均共有n項
    但我希望寫二項您看的懂
    剩下的我再自行想像
    麻煩您起個頭
    拜託了
  • les
  • 抱歉
    其相關的式子我覺得應該說請清楚點
    所以下為補充
    min: Var(rp) = w1^2 * var1 + w2^2 * var2 + 2 * w1 * w2 * sig1 * sig2
    E(rp) = w1 * E(r1) + w2 * E(r2)

    其中 w1+w2 = 1
    var1 = sig1^2 : sig1為E(r1)的標準差
    var2 = sig2^2 : sig2為E(r2)的標準差

    其實變數與權重的部分均共有n項
    但我希望寫二項您看的懂
    剩下的我再自行想像
    麻煩您起個頭
    拜託了
  • 心冷熱情熄
  • 規劃求解...? 依據你的命題來說,常用名詞應該是非線性規劃吧...
    而你列的式子常用在多目標規劃的標準例題,用來說明多目標間各標的互相取捨關係,以利轉換為非線性規劃求最佳解,這種標準命題的解法與圖形在教科書上及相關文獻都有,通常教科書上的範例都只有兩個變數,在你的命題中,var1, var2 看起來是變數,w1, w2 看起來是多目標規劃用的權值,而你有 n 項,該不是想建立通用式吧?
    你的問題已經是數學規劃裡面比較深入的部份,通常是碩士論文級在基礎理論闡述中的命題,並將此命題延伸到個人的論文應用中,程式化的部份通常是 case by case ,目前數學發展尚難建立通式,如果是大學畢業專題,最好換個題目。
    因為這後面跟著的是非線性規劃命題要用哪種最佳解搜尋法,光是這個搜尋法就多的亂七八糟,現在正熱門的遺傳演算法也算是這種,傳統最有名的就是最陡坡降法,還牽扯一維或多維搜尋法,向量或全域搜尋法,而且這個牽扯完整命題才能選擇較合適的最佳解搜尋法,可是從你簡單列出來的部份,看不出來你的數學模型你自己是否足夠了解。
    從命題跟實際可能的運用,在理學院或管理學院,大概是碩士論文級的,在工學院看科系,可能是碩士或博士級的論文,若是碩士論文以上的等級,最好跟你指導教授或學長討論,後面隱藏的非線性規劃最佳解搜尋法是大問題,而這類實際問題求解,一般研究案都至少百萬元級以上,還要做敏感度分析等來調校模式,若是作業級的程式,通常書本或老師會先指定最佳解搜尋法,反而多目標規劃這邊外包裝沒啥特別的,或是你已購買現有的 package 來做後面的最佳解搜尋法?
    你會用規劃求解這四個字,看起來應該對這項學門瞭解不夠深入,最好換個題目。因為這部份是整個數學規劃 (包含線性規劃、動態規劃、整數規劃、網路規劃、非線性規劃等) 可以組成的多目標規劃,要原文書跟相關原文期刊看到一定程度下,才能找出比較合適的方法解決,不然就是改用只比窮舉法稍微好一點的遺傳演算法,反正國內很多學者搞不清楚,還以為你的東西很好... 遺傳演算法實在是沒有辦法中的辦法,但是不須了解整個數學規劃,搞到最後被國內學者當作是無上法寶,所以拿來騙一些外行的審查委員,是很好的方法,才有可能在短時間畢業,否則光是數學規劃理論,起碼搞掉你一年。
  • Shawn
  • 大家共同討論一下.:)
    看起來你是有n個random variables (是以你要用mean and variance來求解)
    X=sigma(Wn*Xn)...而你想要minimize X的variance under the condition
    sigma(Wn)=1...
    一般來說,no constraint optimization有幾種比較常用的方法..大致分為兩類
    1.Taylor expansion based:
    f(x+p)=f(x)+p'*g(x)+p'*H(x)*p+O(x^3)
    where g(x): frist derivative of f(x), H(x): second derivative of f(x)
    Here, I do not limit x has only one element (p' means p transpose)
    based on this expansion, several methods can be derived
    (1) Steepest Descent: (我猜可能是心冷說的最陡坡降法)
    只求g(x),用p=-g(x)..永遠用x-p*stepsize猜下一個x
    (2) Newton Method:
    求g(x) and H(x), 用H(x)*p=-g(x)裡面求出的p作為下一個x的系數
      因為解這個需要用inv(H(x))..所以有很多有的沒的方法就可以用來簡化complexity..像quasi-newton..truncated newton都是很好用的方法
      Newton Method will always outperform steepest decent method..
    but with a high computational cost and complexity
    裡面我忽略了找stepsize的方法.:P 因為很難打字..不過你去查應該都查得到.:))
    2. 用亂猜的...(這樣說有點籠統..不過是真的.:P) 又叫作no-derivative method
    (1) Finite difference: 用現有的點來猜微分值 x'=x(n+1)-x(n)
    (2) Nelder-Meade algorithm: Matlab裡有...就是用幾個點比較.有更小的值就繼續一直作...(就像用猜的.:P,不過它是有用一些幾何在裡面)
    (3) Pattern Search: 這個我不太熟,好像是要找positive basis....
    請用過的人指教.^_^..
    剛我說的都是no constraint下可直接用的...加了constraint之後...
    depends on your constraints, 可以選用barrier method (good for inequality constraint),或是feasible direction formulation(good for equality)..像你的case..你可以直接用feasible direction by assuming
    X=x0+k*[1,-1]'...或是像lagragian...
    套用之後再以前面的方法求optimum....都很好用喔!:)
    建議你可以看看matlab的optimization toolbox,裡面有蠻多介紹的喔.:)
  • les
  • 忑珂
    覜珴媼弇攬衭腔膘荽
    扂慮??艘
    菴媼?陑濮腔隙
    淩岆驚督
    嶒稛岆扂腔幅尪?恅笢腔郔笭珨諔
    ?扂掛砑挲赻撩婓讀讀艘
    祥綎政婓擒褩?嶲怮猧賸
    扂淩腔竭褣
    眕狟扂都岆瞳蚚斕垀枑腔價秪栳呾腔源楊賤
    筍岆?揃家褫眕
    嗣黺?
    惝?宒赽衱?腔竭墿珨揹
    奧嶒悷啃堤
    珨硉蚰祥善螹髯
    秪扂腔翋最宒岆芵綎RATS最宒蚴?腔
    祥眭耋蠟岆瘁衄?綎(耷濮嬡腔?闚,祥岆扂砑蚚腔,硌諒忨?隅腔,漲扂珩賸淕淕圉爛,符類堤珨螹髯,奧坻怮疆,跦掛嶲奪扂侚魂)
    祫黺?墿些
    扂?菴珨竭褫
    衱腕?奧赻撩衱祥宒燴馱掖劓
    婓蚴?腔綎最,珩褫??RATS
    垀眕珨硉婓勛寯嬡輊
    疑祥眢?堤湮窒煦翋最宒
    筍稛賤稛窒煦詻侚
    祫黺蠟枑?觳醴腔膘荽
    扂飲眒?詻善稛窒煦賸飲祥眭耋睡隙斕賸
    祫黺斕垀枑腔淕??? (婦漪?俶??﹜??﹜淕??﹜鋒繚??﹜準?俶??脹) 腕窒煦
    ?扂壺賸鋒繚??婓?砑
    坻飲衄蚚奻
    奧狟醱岆扂?芵綎VBA?傖腔
    ?岆衄悷船
    KEY奻
    洷咡鷂蠟?僕肮??艘艘岆瘁衄賤腔源宒
    祫黺Shawn腔枑荽扂珂梑梑艘蠟垀枑腔matlab髡夔
    螜?岆螹髯
    筍淩腔竭珴珴斕?腔隙

    Function PortfolioWeights(expret, rf, retvec, vcvmat)
    '  Returns the Portfolio Weights (for a given value of expected return)
    '  (special cases : MinVar portfolio if expret = -1, Tangent portfolio if rf > 0)
    Dim a, b, c, d, gi, hi
      Dim i As Integer, n As Integer
      Dim uvec1() As Variant, wtsvec() As Variant
      Dim uvec As Variant, vinvmat As Variant
      Dim l As Variant, m As Variant
      n = Application.Count(retvec)
      ReDim uvec1(n)
      ReDim wtsvec(n)
      For i = 1 To n
        uvec1(i) = 1
      Next i
      uvec = Application.Transpose(uvec1)
      If retvec.Columns.Count > retvec.Rows.Count Then retvec = Application.Transpose(retvec)
      vinvmat = Application.MInverse(vcvmat)
      l = Application.MMult(vinvmat, retvec)
      m = Application.MMult(vinvmat, uvec)
    '  the following three rows have scalar answers derived from array functions
    a = Application.Sum(Application.MMult(Application.Transpose(uvec), l))
      b = Application.Sum(Application.MMult(Application.Transpose(retvec), l))
      c = Application.Sum(Application.MMult(Application.Transpose(uvec), m))
      d = b * c - a ^ 2
    '  minimum variance portfolio
    If expret = -1 Then expret = a / c
    '  optimal risky portfolio / tangent portfolio
    If rf > 0 Then expret = (a / c) - (d / c ^ 2) / (rf - a / c)
      For i = 1 To n
        gi = b * m(i, 1) - a * l(i, 1)
        hi = c * l(i, 1) - a * m(i, 1)
        wtsvec(i) = (gi + hi * expret) / d
      Next i
      PortfolioWeights = wtsvec
    End Function

    狟醱岆扂腔?觳宒赽
    淩洷咡斕?夔艘腔雅扂婓鎊
    Using Algebra to reproduce Unconstrained Frontier Portfolios

    Asset Data Exp Ret Std Dev Portfolio Weights
    TBills 0.6% 4.3% TBills 33.3%
    Bonds 2.1% 10.1% Bonds 33.3%
    Shares 9.0% 20.8% Shares 33.3%

    Correlation Matrix TBills Bonds Shares via fn
    TBills 1.00 0.63 0.09 Exp Ret 3.90% 3.90%
    Bonds 0.63 1.00 0.23 Variance 0.0080
    Shares 0.09 0.23 1.00 Std Dev 8.95% 8.95%

    VCV matrix TBills Bonds Shares VCV inverse
    TBills 0.0018 0.0027 0.0008 901.51 -246.92 10.80
    Bonds 0.0027 0.0102 0.0048 -246.92 171.13 -14.52
    Shares 0.0008 0.0048 0.0433 10.80 -14.52 24.53


    Finding weights, g and h, to generate points on the frontier

    uvec l m g h g+h
    A 3.97
    1 1.20 665.40 B 0.20 124.0% -1851.9% -1727.9%
    1 0.81 -90.30 C 595.91 -20.5% 805.2% 784.6%
    1 1.97 20.82 D 104.15 -3.5% 1046.7% 1043.2%

    Exp Ret 0.0% 100.0% 100.0%
    Std Dev 4.4% 237.6%

    Generating Frontier Portfolios, using g and h

    Target expected return 7.0%

    weights via fn
    TBills -5.6% -5.6% Exp Ret 7.0%
    Bonds 35.8% 35.8% Std Dev 15.7%
    Shares 69.8% 69.8%


    Finding minimum variance portfolio

    Min variance expected return 0.67%

    weights via fn
    TBills 111.7% 111.7% Exp Ret 0.67%
    Bonds -15.2% -15.2% Std Dev 4.10%
    Shares 3.5% 3.5%

    Given existence of risk-free asset, finding optimal risky portfolio
    Risk-free return 0.50%
    Optimal risky expected return 18.32%
    weights via fn
    TBills -215.2% -215.2% Exp Ret 18.32%
    Bonds 127.0% 127.0% Std Dev 42.42%
    Shares 188.2% 188.2%
  • les
  • 天阿 怎麼辯亂碼
    二位大哥對不起喔
    其大概內容就是感謝您們的回應
    與我的問題點
    希望你們能看的懂
  • Shawn
  • sorry..我實在不知道你要弄啥耶.:)
    我看得出來你要作財務計算
    可否請你將你要optimize的equation 寫出,還有你的constraints...
    如果跟上次相同的話...
    constraints: w1+w2+...wn=1
    Equation: Var(w1*x1+w2*x2+....wn*xn)
    -------------------------------------
    你可以用Newton method 加上lagragian 去作...
    當然,所有的numerical method 都不能guarantee會找到global optimum
    但是local optimum 一定沒問題....
  • les
  • 感謝shawn的回應
    您真的好利害
    您提的部分其實我都有試過
    不過可能我有些理論觀念沒搞懂
    所以我現在在補k
    其實我做的有關於投資組合的題目
    大概意識事透過不同風險預測模型的比較
    透過不同資產所預測的風險
    建立投資組合
    min: Var(rp) = w1^2 * var1 + w2^2 * var2 + 2 * w1 * w2 * sig1 * sig2











  • les
  • 感謝shawn的回應
    您與心冷真的都好利害
    您提的部分其實我都有試過
    不過可能我有些理論觀念沒搞懂
    所以我現在在補k
    其實我做的有關於投資組合的題目
    大概意識事透過不同風險預測模型的比較
    透過不同資產所預測的風險
    建立投資組合
    min: Var(rp) = w1^2 * var1 + w2^2 * var2 + 2 * w1 * w2 * sig1 * sig2
    E(rp) = w1 * E(r1) + w2 * E(r2)
    constraints: w1+w2+...wn=1
    w代表的是我的權重
    var1式我的第一項資產的報酬率變異數
    var2式我的第二項資產的報酬率變異數
    var1 = sig1^2
    var2 = sig2^2
    sig 式代表個別資產的標準差
    世子的最後那串是我的共變異數
    其大概是亦是透過給定投資組合風險下,就是投資組合得風險固定下
    求得投資組合的報酬率最高的那一點
    其實這個部分我已經使用另外一套軟體跑出來了
    不過那套軟體的缺點並不方便建立通式
    而且執行中消耗記憶體超級強
    我為了執行這個論文題目
    我弄了四十幾台研究室電腦
    連續跑了快二個星期採跑完
    但是我現在論文有些部分需要修改
    一想到要重跑我就有點放棄
    想說重新寫一個程式可以通用的
    而且不用那麼麻煩的
    不知道二位有沒有聽過RATS與EVIEWS
    其實二套適用於跑統計
    我盡然笨笨的拿來跑計算
    而且我盡然賭氣拿RATS跟教授挑戰
    已經不知道要說啥了
    真的有點小後悔
    但是畢竟有些地方真的也很好用
    其實教授已經不太知道我怎弄得拉
    所以我也沒辦法解決那個煩人的跑程式的時間
    所以我才想與大家討論看看
    我把我論文中的主程式波上來
    大家來看看是否我的問題還有沒解決的
    不然怎會跑那麼久













  • les
  • 首先先說明
    因為我看了很多之前的文章
    大多上來求作業
    但是我並沒有那個意識
    所以我將我的主程式拿出來與各位討論
    希望大家別誤會

    *******************************************************************
       計算PACIFIC的四種比較方法 (歷史) 以美元計價
    給定風險下投資組合之歷史資料報酬率最高次數的結果 ***
    ************************************************************
    all 162
    open data return.wks
    data(for=wks,org=col) / a_r n_r j_r h_r s_r
    ********************************** 定義變異數資料總數           **
    *compute n_var  = 127
    ********************************** 定義歷史資料大迴圈(模擬次數)      **
    compute n_draw = 1
    **********************************   定義方法種類              **
    compute methods = 4
    ********************************** 定義投資組合中的國家數目        **
    compute countries = 5
    compute nations = countries
    ********************************** 定義報酬率資料總數         **
    compute n_obs = 162
    *compute n_obs = 10
    ********************************** 定義前置資料總數         **
    compute n_win =  3
    compute n_var  = n_obs - n_win - 1
    ********************************** 定義指數平滑lambda值         **
    compute lambda03 = 0.03
    ********************************** 定義各方法最大及最小值次數      **
    compute all_max_count = 0
    compute f03_max_count = 0
    compute e03_max_count = 0
    compute gar_max_count = 0
    compute all_min_count = 0
    compute f03_min_count = 0
    compute e03_min_count = 0
    compute gar_min_count = 0
    ********************************** 定義計算報酬率相關矩陣         **
    dec rec all_sigma(1,nations)
    dec rec f03_sigma(1,nations)
    dec rec e03_sigma(1,nations)
    dec rec gar_sigma(1,nations)
    dec rec all_d(1,countries)
    dec rec f03_d(1,countries)
    dec rec e03_d(1,countries)
    dec rec gar_d(1,countries)
    dec rec all_w(1,countries)
    dec rec f03_w(1,countries)
    dec rec e03_w(1,countries)
    dec rec gar_w(1,countries)
    dec rec vm_mat(n_var,methods)
    dec rec rp_mat(n_var,methods)
    dec rec r_real_mat(countries,n_var)
    ********************************** 定義投資組合相關矩陣         **
    dec rec all_amat(countries,countries)
    dec rec f03_amat(countries,countries)
    dec rec e03_amat(countries,countries)
    dec rec gar_amat(countries,countries)
    dec rec all_c_mat(1,countries)
    dec rec f03_c_mat(1,countries)
    dec rec e03_c_mat(1,countries)
    dec rec gar_c_mat(1,countries)
    dec rec all_cvmat(countries,countries)
    dec rec f03_cvmat(countries,countries)
    dec rec e03_cvmat(countries,countries)
    dec rec gar_cvmat(countries,countries)
    dec rec all_c(countries+1,countries+1)
    dec rec f03_c(countries+1,countries+1)
    dec rec e03_c(countries+1,countries+1)
    dec rec gar_c(countries+1,countries+1)
    ********************************** 定義 vm 相關矩陣         **
    dec rec vm_var(methods,n_var)
    ********************************** 定義變異數相關矩陣         **
    dec rec a_all_var(n_draw,n_var)
    dec rec a_f03_var(n_draw,n_var)
    dec rec a_e03_var(n_draw,n_var)
    dec rec a_gar_var(n_draw,n_var)
    dec rec n_all_var(n_draw,n_var)
    dec rec n_f03_var(n_draw,n_var)
    dec rec n_e03_var(n_draw,n_var)
    dec rec n_gar_var(n_draw,n_var)
    dec rec j_all_var(n_draw,n_var)
    dec rec j_f03_var(n_draw,n_var)
    dec rec j_e03_var(n_draw,n_var)
    dec rec j_gar_var(n_draw,n_var)
    dec rec h_all_var(n_draw,n_var)
    dec rec h_f03_var(n_draw,n_var)
    dec rec h_e03_var(n_draw,n_var)
    dec rec h_gar_var(n_draw,n_var)
    dec rec s_all_var(n_draw,n_var)
    dec rec s_f03_var(n_draw,n_var)
    dec rec s_e03_var(n_draw,n_var)
    dec rec s_gar_var(n_draw,n_var)
    ********************************** 指定gar序列              **
    declare series u
    declare series v
    ********************************** 進行模擬迴圈         **
    * compute n_simu = 0
    display @11 "*** 開始運算歷史資料 ***"
    * while n_simu < n_draw {
    compute n_simu = 1
    **********************************   計算PACIFIC地區的預測變異數      **
    * compute n_simu = n_simu +1
    * display "n_simu" n_simu
    **********************************   計算澳洲預測變數
    **********************************   GARCH (A)
    * loop
    * compute gar_temp = -1.0
    * until gar_temp < 0
    * boot entries / 1 n_obs
    * set a_r = aord(entries(t))
    do j = 1,n_var
    set a_gar_r 1 n_win+j = a_r
    compute start=1
    set u = 0.0; set v = 0.0;
    nonlin b0 b1 a0 a1 a2
    frml regresid = a_gar_r-b0-b1*a_gar_r{1}
    frml garchvar = a0+a1*u{1}**2+a2*v{1}
    frml garchlogl = v(t)=garchvar(t),u(t)=regresid(t),-.5*(log(v(t)) + (u(t)**2/v(t)))
    linreg(noprint) a_gar_r / u
    # constant a_gar_r{1}
    nlpar(subiter=1000)
    compute b0=%beta(1),b1=%beta(2)
    compute a0=%seesq,a1=0.1,a2=-.2
    set v =%seesq
    maximize(method=simplex,recursive,iter=6,noprint) garchlogl start+1 n_win+j
    maximize(method=bhhh,recursive,iterations=10000,noprint) garchlogl start+1 n_win+j
    compute a_gar_var(n_simu,j) = sqrt(abs(a0+a1*u(j+n_win)**2+a2*v(j+n_win)))
    end do j
    ***************** 以下為測試時要顯示的
    * write a_gar_var
    * display "a0" a0
    * display "a1" a1
    * display "u" j+n_win u(j+n_win)
    * display "v" j+n_win v(j+n_win)
    **************** 以下為模擬用
    * compute gar_temp = a0+a1*u(j+n_win)**2+a2*v(j+n_win)
    * display "a_garchvar" n_simu j gar_temp
    * if gar_temp < 0 {
    **compute n_simu = n_simu - 1
    * break
    * }
    * else
    * compute a_gar_var(n_simu,j) = sqrt(a0+a1*u(j+n_win)**2+a2*v(j+n_win))
    * end do j
    * if gar_temp > 0
    * break
    * end loop
    ********************************** 計算紐西蘭預測變異數
    **********************************   GARCH (N)
    * loop
    * compute gar_temp = -1.0
    * until gar_temp < 0
    * boot entries / 1 n_obs
    * set n_r = nzse(entries(t))
    do j = 1,n_var
    set n_gar_r 1 n_win+j = n_r
    compute start=1
    set u = 0.0; set v = 0.0;
    nonlin b0 b1 a0 a1 a2
    frml regresid = n_gar_r-b0-b1*n_gar_r{1}
    frml garchvar = a0+a1*u{1}**2+a2*v{1}
    frml garchlogl = v(t)=garchvar(t),u(t)=regresid(t),-.5*(log(v(t)) + (u(t)**2/v(t)))
    linreg(noprint) n_gar_r / u
    # constant n_gar_r{1}
    nlpar(subiter=1000)
    compute b0=%beta(1),b1=%beta(2)
    compute a0=%seesq,a1=0.1,a2=-.2
    set v =%seesq
    maximize(method=simplex,recursive,iter=6,noprint) garchlogl start+1 n_win+j
    maximize(method=bhhh,recursive,iterations=10000,noprint) garchlogl start+1 n_win+j
    compute n_gar_var(n_simu,j) = sqrt(abs(a0+a1*u(j+n_win)**2+a2*v(j+n_win)))
    end do j
    ***************** 以下為測試時要顯示的
    * write a_gar_var
    * display "a0" a0
    * display "a1" a1
    * display "u" j+n_win u(j+n_win)
    * display "v" j+n_win v(j+n_win)
    **************** 以下為模擬用
    * compute gar_temp = a0+a1*u(j+n_win)**2+a2*v(j+n_win)
    * display "n_garchvar" n_simu j gar_temp
    * if gar_temp < 0 {
    **compute n_simu = n_simu - 1
    * break
    * }
    * else
    * compute n_gar_var(n_simu,j) = sqrt(a0+a1*u(j+n_win)**2+a2*v(j+n_win))
    * end do j
    * if gar_temp > 0
    * break
    * end loop
    ********************************** 計算日本預測變異數
    **********************************   GARCH (J)
    * loop
    * compute gar_temp = -1.0
    * until gar_temp < 0
    * boot entries / 1 n_obs
    * set j_r = jp(entries(t))
    do j = 1,n_var
    set j_gar_r 1 n_win+j = j_r
    compute start=1
    set u = 0.0; set v = 0.0;
    nonlin b0 b1 a0 a1 a2
    frml regresid = j_gar_r-b0-b1*j_gar_r{1}
    frml garchvar = a0+a1*u{1}**2+a2*v{1}
    frml garchlogl = v(t)=garchvar(t),u(t)=regresid(t),-.5*(log(v(t)) + (u(t)**2/v(t)))
    linreg(noprint) j_gar_r / u
    # constant j_gar_r{1}
    nlpar(subiter=1000)
    compute b0=%beta(1),b1=%beta(2)
    compute a0=%seesq,a1=0.1,a2=-.2
    set v =%seesq
    maximize(method=simplex,recursive,iter=6,noprint) garchlogl start+1 n_win+j
    maximize(method=bhhh,recursive,iterations=10000,noprint) garchlogl start+1 n_win+j
    compute j_gar_var(n_simu,j) = sqrt(abs(a0+a1*u(j+n_win)**2+a2*v(j+n_win)))
    end do j
    ***************** 以下為測試時要顯示的
    * write a_gar_var
    * display "a0" a0
    * display "a1" a1
    * display "u" j+n_win u(j+n_win)
    * display "v" j+n_win v(j+n_win)
    **************** 以下為模擬用
    * compute gar_temp = a0+a1*u(j+n_win)**2+a2*v(j+n_win)
    **display "j_garchvar" n_simu j gar_temp
    * if gar_temp < 0 {
    **compute n_simu = n_simu - 1
    * break
    * }
    * else
    * compute j_gar_var(n_simu,j) = sqrt(a0+a1*u(j+n_win)**2+a2*v(j+n_win))
    * end do j
    * if gar_temp > 0
    * break
    * end loop
    ********************************** 計算香港預測變數
    **********************************   GARCH (HSI)
    do j = 1,n_var
    set h_gar_r 1 n_win+j = h_r
    compute start=1
    set u = 0.0; set v = 0.0;
    nonlin b0 b1 a0 a1 a2
    frml regresid = h_gar_r-b0-b1*h_gar_r{1}
    frml garchvar = a0+a1*u{1}**2+a2*v{1}
    frml garchlogl = v(t)=garchvar(t),u(t)=regresid(t),-.5*(log(v(t)) + (u(t)**2/v(t)))
    linreg(noprint) h_gar_r / u
    # constant h_gar_r{1}
    nlpar(subiter=1000)
    compute b0=%beta(1),b1=%beta(2)
    compute a0=%seesq,a1=0.1,a2=-.2
    set v =%seesq
    maximize(method=simplex,recursive,iter=6,noprint) garchlogl start+1 n_win+j
    maximize(method=bhhh,recursive,iterations=10000,noprint) garchlogl start+1 n_win+j
    compute h_gar_var(n_simu,j) = sqrt(abs(a0+a1*u(j+n_win)**2+a2*v(j+n_win)))
    end do j
    ********************************** 計算新加坡預測變異數
    **********************************   GARCH (STI)
    do j = 1,n_var
    set s_gar_r 1 n_win+j = s_r
    compute start=1
    set u = 0.0; set v = 0.0;
    nonlin b0 b1 a0 a1 a2
    frml regresid = s_gar_r-b0-b1*s_gar_r{1}
    frml garchvar = a0+a1*u{1}**2+a2*v{1}
    frml garchlogl = v(t)=garchvar(t),u(t)=regresid(t),-.5*(log(v(t)) + (u(t)**2/v(t)))
    linreg(noprint) s_gar_r / u
    # constant s_gar_r{1}
    nlpar(subiter=1000)
    compute b0=%beta(1),b1=%beta(2)
    compute a0=%seesq,a1=0.1,a2=-.2
    set v =%seesq
    maximize(method=simplex,recursive,iter=6,noprint) garchlogl start+1 n_win+j
    maximize(method=bhhh,recursive,iterations=10000,noprint) garchlogl start+1 n_win+j
    compute s_gar_var(n_simu,j) = sqrt(abs(a0+a1*u(j+n_win)**2+a2*v(j+n_win)))
    end do j
    **********************************   (澳洲)
    ********************************** All Data
    do j = 1,n_var
    stat(noprint) a_r 2 n_win+j
    compute a_all_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Fix-03
    do j = 1,n_var
    stat(noprint) a_r j+1 n_win+j
    compute a_f03_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Ex-03
    *set a_e03_var 1 n_var = 0.0
    do j=1,n_var
    set a_e03_r 1 n_win+j = a_r
    acc a_e03_r / a_sum
    set a_e03_r_t = a_sum / t
    set a_e03_error2 = (a_e03_r - a_e03_r_t)**2
    esmooth(alpha=0.97,smoothed=a_e03_temp_var) a_e03_error2 2 n_win+j
    compute a_e03_var(n_simu,j) = sqrt((1-lambda03)*a_e03_error2(j+n_win)+lambda03*a_e03_temp_var(j+n_win))
    end do j
    ********************************** (紐西蘭)
    ********************************** All Data
    do j = 1,n_var
    stat(noprint) n_r 2 n_win+j
    compute n_all_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Fix-03
    do j = 1,n_var
    stat(noprint) n_r j+1 n_win+j
    compute n_f03_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Ex-03
    *set n_e03_var 1 n_var = 0.0
    do j=1,n_var
    set n_e03_r 1 n_win+j = n_r
    acc n_e03_r / n_sum
    set n_e03_r_t = n_sum / t
    set n_e03_error2 = (n_e03_r - n_e03_r_t)**2
    esmooth(alpha=0.97,smoothed=n_e03_temp_var) n_e03_error2 2 n_win+j
    compute n_e03_var(n_simu,j) = sqrt((1-lambda03)*n_e03_error2(j+n_win)+lambda03*n_e03_temp_var(j+n_win))
    end do j
    ********************************** (日本)
    ********************************** All Data
    do j = 1,n_var
    stat(noprint) j_r 2 n_win+j
    compute j_all_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Fix-03
    do j = 1,n_var
    stat(noprint) j_r j+1 n_win+j
    compute j_f03_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Ex-03
    *set j_e03_var 1 n_var = 0.0
    do j=1,n_var
    set j_e03_r 1 n_win+j = j_r
    acc j_e03_r / j_sum
    set j_e03_r_t = j_sum / t
    set j_e03_error2 = (j_e03_r - j_e03_r_t)**2
    esmooth(alpha=0.97,smoothed=j_e03_temp_var) j_e03_error2 2 n_win+j
    compute j_e03_var(n_simu,j) = sqrt((1-lambda03)*j_e03_error2(j+n_win)+lambda03*j_e03_temp_var(j+n_win))
    end do j
    **********************************   (香港)
    ********************************** All Data
    do j = 1,n_var
    stat(noprint) h_r 2 n_win+j
    compute h_all_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Fix-03
    do j = 1,n_var
    stat(noprint) h_r j+1 n_win+j
    compute h_f03_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Ex-03
    *set h_e03_var 1 n_var = 0.0
    do j=1,n_var
    set h_e03_r 1 n_win+j = h_r
    acc h_e03_r / h_sum
    set h_e03_r_t = h_sum / t
    set h_e03_error2 = (h_e03_r - h_e03_r_t)**2
    esmooth(alpha=0.97,smoothed=h_e03_temp_var) h_e03_error2 2 n_win+j
    compute h_e03_var(n_simu,j) = sqrt((1-lambda03)*h_e03_error2(j+n_win)+lambda03*h_e03_temp_var(j+n_win))
    end do j
    **********************************   (新加坡)
    ********************************** All Data
    do j = 1,n_var
    stat(noprint) s_r 2 n_win+j
    compute s_all_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Fix-03
    do j = 1,n_var
    stat(noprint) s_r j+1 n_win+j
    compute s_f03_var(n_simu,j) = sqrt(%variance)
    end do j
    ********************************** Ex-03
    *set s_e03_var 1 n_var = 0.0
    do j=1,n_var
    set s_e03_r 1 n_win+j = s_r
    acc s_e03_r / s_sum
    set s_e03_r_t = s_sum / t
    set s_e03_error2 = (s_e03_r - s_e03_r_t)**2
    esmooth(alpha=0.97,smoothed=s_e03_temp_var) s_e03_error2 2 n_win+j
    compute s_e03_var(n_simu,j) = sqrt((1-lambda03)*s_e03_error2(j+n_win)+lambda03*s_e03_temp_var(j+n_win))
    end do j
    ********************************** 給定預期風險下投資組合之最大報酬率   **
    ********************************** 定義資料起始點         **
    compute start = 1
    compute end  = n_obs
    ********************************** 計算各國報酬率平均值         **
    accumulate a_r / sum_a_r
    compute a_av_r = sum_a_r(end)/end
    accumulate n_r / sum_n_r
    compute n_av_r = sum_n_r(end)/end
    accumulate j_r / sum_j_r
    compute j_av_r = sum_j_r(end)/end
    accumulate h_r / sum_h_r
    compute h_av_r = sum_h_r(end)/end
    accumulate s_r / sum_s_r
    compute s_av_r = sum_s_r(end)/end
    ********************************** 建立平均報酬率矩陣(1×5)        **
    compute av_r = !!a_av_r,n_av_r,j_av_r,h_av_r,s_av_r!!
    ********************************** 建立相關係數矩陣          **
    cmom(corr,noprint)
    # a_r n_r j_r h_r s_r
    ********************************** 建立變異數矩陣             **
    set a_all_temp 1 n_var = a_all_var(n_simu,t)
    set n_all_temp 1 n_var = n_all_var(n_simu,t)
    set j_all_temp 1 n_var = j_all_var(n_simu,t)
    set h_all_temp 1 n_var = h_all_var(n_simu,t)
    set s_all_temp 1 n_var = s_all_var(n_simu,t)
    make all_mat_var start end
    # a_all_temp n_all_temp j_all_temp h_all_temp s_all_temp
    set a_f03_temp 1 n_var = a_f03_var(n_simu,t)
    set n_f03_temp 1 n_var = n_f03_var(n_simu,t)
    set j_f03_temp 1 n_var = j_f03_var(n_simu,t)
    set h_f03_temp 1 n_var = h_f03_var(n_simu,t)
    set s_f03_temp 1 n_var = s_f03_var(n_simu,t)
    make f03_mat_var start end
    # a_f03_temp n_f03_temp j_f03_temp h_f03_temp s_f03_temp
    set a_e03_temp 1 n_var = a_e03_var(n_simu,t)
    set n_e03_temp 1 n_var = n_e03_var(n_simu,t)
    set j_e03_temp 1 n_var = j_e03_var(n_simu,t)
    set h_e03_temp 1 n_var = h_e03_var(n_simu,t)
    set s_e03_temp 1 n_var = s_e03_var(n_simu,t)
    make e03_mat_var start end
    # a_e03_temp n_e03_temp j_e03_temp h_e03_temp s_e03_temp
    set a_gar_temp 1 n_var = a_gar_var(n_simu,t)
    set n_gar_temp 1 n_var = n_gar_var(n_simu,t)
    set j_gar_temp 1 n_var = j_gar_var(n_simu,t)
    set h_gar_temp 1 n_var = h_gar_var(n_simu,t)
    set s_gar_temp 1 n_var = s_gar_var(n_simu,t)
    make gar_mat_var start end
    # a_gar_temp n_gar_temp j_gar_temp h_gar_temp s_gar_temp
    ********************************** 抓各種方法的5個國家各月變異數     **
    set all_riskp start n_var = 0.0
    set f03_riskp start n_var = 0.0
    set e03_riskp start n_var = 0.0
    set gar_riskp start n_var = 0.0
    *********************************    抓C矩陣的共變異數的部份         **
    ********************************* 定義實際報酬率矩陣
    set real_a_r n_win+2 n_obs = a_r
    set real_n_r n_win+2 n_obs = n_r
    set real_j_r n_win+2 n_obs = j_r
    set real_h_r n_win+2 n_obs = h_r
    set real_s_r n_win+2 n_obs = s_r

    make r_temp_mat n_win+2 n_obs
    # real_a_r real_n_r real_j_r real_h_r real_s_r
    compute r_real_mat = tr(r_temp_mat)
    do months=1,n_var
    do nations=1,countries
    compute all_sigma(1,nations) = all_mat_var(months,nations)
    compute f03_sigma(1,nations) = f03_mat_var(months,nations)
    compute e03_sigma(1,nations) = e03_mat_var(months,nations)
    compute gar_sigma(1,nations) = gar_mat_var(months,nations)
    end do nations
    ********************************** 建立各月C矩陣         **
    ****************** ALL **
    do i=1,countries+1
    do j=1,countries+1
    if i<(countries+1).and.j<(countries+1)
    compute all_c(i,j) = %cmom(i,j)*all_sigma(1,i)*all_sigma(1,j)*2
    else if i<>j
    compute all_c(i,j) = 1
    end do j
    end do i
    compute all_c(countries+1,countries+1)=0
    ****************** F03 **
    do i=1,countries+1
    do j=1,countries+1
    if i<(countries+1).and.j<(countries+1)
    compute f03_c(i,j) = %cmom(i,j)*f03_sigma(1,i)*f03_sigma(1,j)*2
    else if i<>j
    compute f03_c(i,j) = 1
    end do j
    end do i
    compute f03_c(countries+1,countries+1)=0
    ****************** E03 **
    do i=1,countries+1
    do j=1,countries+1
    if i<(countries+1).and.j<(countries+1)
    compute e03_c(i,j) = %cmom(i,j)*e03_sigma(1,i)*e03_sigma(1,j)*2
    else if i<>j
    compute e03_c(i,j) = 1
    end do j
    end do i
    compute e03_c(countries+1,countries+1)=0
    ****************** GAR **
    do i=1,countries+1
    do j=1,countries+1
    if i<(countries+1).and.j<(countries+1)
    compute gar_c(i,j) = %cmom(i,j)*gar_sigma(1,i)*gar_sigma(1,j)*2
    else if i<>j
    compute gar_c(i,j) = 1
    end do j
    end do i
    compute gar_c(countries+1,countries+1)=0
    ********************************** 建立各月C之反矩陣         **
    compute all_invC = inv(all_c)
    compute f03_invC = inv(f03_c)
    compute e03_invC = inv(e03_c)
    compute gar_invC = inv(gar_c)
    **********************************************
    * 建立amat矩陣 ( 為反矩陣中前 5 行 5 列 )***
    **********************************************
    do i=1,countries
    do j=1,countries
    compute all_amat(i,j) = all_invC(i,j)
    compute f03_amat(i,j) = f03_invC(i,j)
    compute e03_amat(i,j) = e03_invC(i,j)
    compute gar_amat(i,j) = gar_invC(i,j)
    end do j
    end do i
    *********************************** 計算d矩陣     **
    compute all_d = av_r * all_amat
    compute f03_d = av_r * f03_amat
    compute e03_d = av_r * e03_amat
    compute gar_d = av_r * gar_amat
    *write all_d
    ********************************** 計算各月KE值              **
    compute tr_r = tr(av_r)
    compute all_ke = %scalar(av_r * all_amat * tr_r)
    compute f03_ke = %scalar(av_r * f03_amat * tr_r)
    compute e03_ke = %scalar(av_r * e03_amat * tr_r)
    compute gar_ke = %scalar(av_r * gar_amat * tr_r)
    *display all_ke
    ********************************** 計算各種方法的各月EM值         **
    do j=1,countries
    compute all_c_mat(1,j) = all_invC(countries+1,j)
    compute f03_c_mat(1,j) = f03_invC(countries+1,j)
    compute e03_c_mat(1,j) = e03_invC(countries+1,j)
    compute gar_c_mat(1,j) = gar_invC(countries+1,j)
    end do j
    compute all_EM = %scalar(all_c_mat * tr_r)
    compute f03_EM = %scalar(f03_c_mat * tr_r)
    compute e03_EM = %scalar(e03_c_mat * tr_r)
    compute gar_EM = %scalar(gar_c_mat * tr_r)
    ********************************** 計算各種方法的各月VM值         **
    do i=1,countries
    do j=1,countries
    compute all_cvmat(i,j) = all_c(i,j)
    compute f03_cvmat(i,j) = f03_c(i,j)
    compute e03_cvmat(i,j) = e03_c(i,j)
    compute gar_cvmat(i,j) = gar_c(i,j)
    end do j
    end do i
    ************************** all data ***
    compute  all_Summe1 = 0.0
    compute  all_Summe2 = 0.0
    compute i = 1
    compute j = 1
    compute temp = 0.0
    do i=1,countries
      do j=i,countries
       IF (i=j)
       compute all_Summe1 = all_Summe1 + (all_c_mat(1,i)*all_c_mat(1,j)*all_cvmat(i,j)/2)
      end do j
     end do i

    do i=1,countries
      do j=1,countries
       IF (i<>j)
        compute all_Summe2 = all_Summe2 + (all_c_mat(1,i)*all_c_mat(1,j)*all_cvmat(i,j))
      end do j
     end do i

    compute all_Summe3 = all_Summe2/2
    compute all_Summe = all_Summe1+all_Summe3
    compute all_vm = %scalar(all_Summe)
    ************************** fixed win03 ***
    compute  f03_Summe1 = 0.0
    compute  f03_Summe2 = 0.0
    compute i = 1
    compute j = 1
    compute temp = 0.0
    do i=1,countries
      do j=i,countries
       IF (i=j)
        compute f03_Summe1 = f03_Summe1 + (f03_c_mat(1,i)*f03_c_mat(1,j)*f03_cvmat(i,j)/2)
       end do j
     end do i

    do i=1,countries
      do j=1,countries
       IF (i<>j)
        compute f03_Summe2 = f03_Summe2 + (f03_c_mat(1,i)*f03_c_mat(1,j)*f03_cvmat(i,j))
      end do j
     end do i

    compute f03_Summe3 = f03_Summe2/2
    compute f03_Summe = f03_Summe1+f03_Summe3
    compute f03_vm = %scalar(f03_Summe)
    ************************** es03 ***
    compute  e03_Summe1 = 0.0
    compute  e03_Summe2 = 0.0
    compute i = 1
    compute j = 1
    compute temp = 0.0
    do i=1,countries
      do j=i,countries
       IF (i=j)
        compute e03_Summe1 = e03_Summe1 + (e03_c_mat(1,i)*e03_c_mat(1,j)*e03_cvmat(i,j)/2)
       end do j
    end do i

    do i=1,countries
      do j=1,countries
       IF (i<>j)
        compute e03_Summe2 = e03_Summe2 + (e03_c_mat(1,i)*e03_c_mat(1,j)*e03_cvmat(i,j))
      end do j
    end do i

    compute e03_Summe3 = e03_Summe2/2
    compute e03_Summe = e03_Summe1+e03_Summe3
    compute e03_vm = %scalar(e03_Summe)
    ************************** garch ***
    compute  gar_Summe1 = 0.0
    compute  gar_Summe2 = 0.0
    compute i = 1
    compute j = 1
    compute temp = 0.0
    do i=1,countries
      do j=i,countries
       IF (i=j)
        compute gar_Summe1 = gar_Summe1 + (gar_c_mat(1,i)*gar_c_mat(1,j)*gar_cvmat(i,j)/2)
       end do j
    end do i

    do i=1,countries
      do j=1,countries
       IF (i<>j)
        compute gar_Summe2 = gar_Summe2 + (gar_c_mat(1,i)*gar_c_mat(1,j)*gar_cvmat(i,j))
      end do j
    end do i

    compute gar_Summe3 = gar_Summe2/2
    compute gar_Summe = gar_Summe1+gar_Summe3
    compute gar_vm = %scalar(gar_Summe)
    ********************************** 比較各月VM 值          **
    compute vm_mat(months,1) = %scalar(all_vm)
    compute vm_mat(months,2) = %scalar(f03_vm)
    compute vm_mat(months,3) = %scalar(e03_vm)
    compute vm_mat(months,4) = %scalar(gar_vm)
    set vm_series 1 methods = 0.0
    do k=1,methods
    compute vm_series(k) = vm_mat(months,k)
    end do k
    ext(noprint) vm_series
    if %maxent == 1
    compute maxVM = all_vm
    else if %maxent == 2
    compute maxVM = f03_vm
    else if %maxent == 3
    compute maxVM = e03_vm
    else
    compute maxVM = gar_vm
    ********************************** 計算各月theta值         **
    *display maxVM
    compute all_theta = %sqrt(2*(maxVM-all_vm)/all_ke)
    compute f03_theta = %sqrt(2*(maxVM-f03_vm)/f03_ke)
    compute e03_theta = %sqrt(2*(maxVM-e03_vm)/e03_ke)
    compute gar_theta = %sqrt(2*(maxVM-gar_vm)/gar_ke)
    *display gar_theta
    *break
    ********************************** 計算各國家權重          **
    do j=1,countries
    compute all_w(1,j) = all_c_mat(1,j) + %scalar(all_theta) * all_d(1,j)
    compute f03_w(1,j) = f03_c_mat(1,j) + %scalar(f03_theta) * f03_d(1,j)
    compute e03_w(1,j) = e03_c_mat(1,j) + %scalar(e03_theta) * e03_d(1,j)
    compute gar_w(1,j) = gar_c_mat(1,j) + %scalar(gar_theta) * gar_d(1,j)
    end do j
    *write all_w
    *write f03_w
    *write e03_w
    *write gar_w
    *write gar_d
    *break
    *display %scalar(gar_theta)*gar_d(1,2)+gar_c_mat(1,2)
    *display gar_theta* gar_d(1,1)
    *display %scalar(gar_theta)*gar_d(1,2)
    *display gar_c_mat(1,2)
    *display %scalar(gar_theta)
    ********************************** 計算各種方法投資組合實際報酬率    **
    compute all_rp = 0.0
    compute f03_rp = 0.0
    compute e03_rp = 0.0
    compute gar_rp = 0.0
    do j = 1,countries
    compute all_rp = all_w(1,j) * r_real_mat(j,months) + all_rp
    compute f03_rp = f03_w(1,j) * r_real_mat(j,months) + f03_rp
    compute e03_rp = e03_w(1,j) * r_real_mat(j,months) + e03_rp
    compute gar_rp = gar_w(1,j) * r_real_mat(j,months) + gar_rp
    end do j
    *display all_w(1,3)
    *display r_real_mat(1,5)
    *display "all_rp" all_rp
    *display "f03_rp" f03_rp
    *display "e03_rp" e03_rp
    *display "gar_rp" gar_rp
    *write r_real_mat
    compute rp_mat(months,1) = %scalar(all_rp)
    compute rp_mat(months,2) = %scalar(f03_rp)
    compute rp_mat(months,3) = %scalar(e03_rp)
    compute rp_mat(months,4) = %scalar(gar_rp)
    ********************************** 累計最高報酬率次數         **
    set rp_series 1 methods = 0.0
    do k=1,methods
    compute rp_series(k) = rp_mat(months,k)
    end do k
    ext(noprint) rp_series
    if %maxent == 1
    compute all_max_count = all_max_count + 1
    else if %maxent == 2
    compute f03_max_count = f03_max_count + 1
    else if %maxent == 3
    compute e03_max_count = e03_max_count + 1
    else
    compute gar_max_count = gar_max_count + 1

    display "             "
    display @2 "***  運算PACIFIC地區的各種方法   ***"
    display @2 "***  計算歷史資料報酬率最高次數的結果  ***"
    display @16 "No_Max"
    display @5 "All" @15 ##### all_max_count
    display @5 "F03" @15 ##### f03_max_count
    display @5 "E03" @15 ##### e03_max_count
    display @5 "Gar" @15 ##### gar_max_count
    display "---------------------------------------"
    end do months
    ***************各種方法實際報酬率呈現
    *write rp_mat

    *print / real_a_r real_n_r real_j_r
    *write r_temp_mat
    *display r_real_mat(3,8)


You haven’t logged in yet, please use guest status to leave message. You can also log in with above service account and leave message

other options